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Abstract 

Thermospheric  density  impacts  satellite  position  and  lifetime  through  atmospheric 
drag.  More  accurate  specification  of  thermospheric  temperature,  a  key  input  to  current 
models  such  as  the  High  Accuracy  Satellite  Drag  Model  (HASDM),  can  decrease  model 
density  errors.  This  thesis  builds  on  Burke  et  al.’s  driven-dissipative  model  (2009)  to 
model  the  arithmetic  mean  temperature,  T\a  ,  defined  by  the  Jacchia,  1977  model  (J77), 
using  the  magnetospheric  electric  field  as  a  driver.  Three  methods  of  treating  the  UV 
contribution  to  Tm  (Ti/2uv)  are  tested.  Two  model  parameters,  the  coupling  and 
relaxation  constants,  are  adjusted  for  38  storms  from  2002  -  2008  to  minimize  modeled 
T 1/2  errors.  Observed  T1/2  values  are  derived  from  densities  and  heights  measured  by  the 
GRACE  satellite.  It  is  found  that  allowing  Ti/2uv  to  vary  produces  the  lowest  errors  for 
27  of  38  storms  in  the  sample  and  27  of  28  storms  with  decreasing  ETV  contributions  over 
the  storm  period.  Treating  Ti/2uv  as  a  constant  produces  the  lowest  errors  for  7  of  10 
storms  with  increasing  UV  contributions.  The  coupling  and  relaxation  constants  were 
found  to  vary  over  the  solar  cycle  and  are  fit  well  as  quadratic  functions  of  V^iova-  By 
using  the  J77  model  to  convert  the  model  T1/2  values  to  density  values,  the  driven- 
dissipative  model  produces  density  errors  slightly  lower  than  HASDM  storm  time  errors. 
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MODELING  THE  THERMOSPHERE  AS  A  DRIVEN-DISSIPATIVE 


THERMODYNAMIC  SYSTEM 


I.  Introduction 


Motivation 

The  thermosphere  is  defined  as  the  neutral  part  of  the  Earth’s  upper  atmosphere 
from  roughly  95  -  1000  km  above  sea  level.  Hundreds  of  Department  of  Defense  and 
other  low-Earth  orbit  satellites  operate  at  these  altitudes.  The  ability  to  accurately 
characterize  the  thermospheric  environment  is  critical  in  an  era  when  the  Department  of 
Defense’s  dependence  on  satellites  for  communications,  intelligence  and  other 
capabilities  has  never  been  higher.  Likewise  as  the  thermosphere  becomes  more  crowded 
with  low-Earth  orbit  satellites  and  space  debris  the  consequences  of  inaccurate  forecasts 
are  becoming  more  significant.  Several  recent  events  illustrate  these  consequences.  The 
destruction  of  the  defunct  Feng  Yun  1C  satellite  by  an  anti-satellite  weapons  test  in 
January,  2007  resulted  in  more  than  2500  new  pieces  of  debris  in  low  earth  orbit  (Burke, 
et  al.,  2009).  The  risk  posed  to  operational  satellites  by  space  debris  was  illustrated  in 
2009  when  the  Iridium  33  satellite  was  destroyed  by  a  collision  with  the  non-operational 
Cosmos  2251  satellite  (Burke,  et  al.,  2010).  There  have  been  several  instances,  such  as  12 
March  and  1  December,  2009,  where  the  risk  of  collision  with  debris  has  forced  the  crew 
of  the  International  Space  Station  to  take  emergency  actions  to  ensure  their  safety 
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(Weimer,  et  al.,  2011).  Improved  characterization  of  the  thermospheric  environment  is 
necessary  to  increase  space  object  tracking  accuracy  and  allow  satellite  operators  and 
manned  spaceflight  missions  to  anticipate  and  avoid  collisions  (Wright,  2007). 


Satellite  Drag 

Variations  in  thermospheric  density  impact  satellite  orbit  trajectories  through 
increased  drag.  The  acceleration  due  to  atmospheric  drag  is  given  by 


where  Asc  and  Msc  are  the  cross-sectional  area  and  mass  of  the  spacecraft,  respectively,  p 
is  the  neutral  mass  density  of  the  atmosphere,  and  V  is  the  spacecraft  velocity  relative  to 
the  neutral  atmosphere.  The  drag  coefficient  CD  depends  on  the  angle  of  flow  to  the 
spacecraft  surface,  the  ratio  of  the  temperatures  of  the  spacecraft  surface  and  the  local 
atmosphere,  and  the  ratio  of  the  mean  mass  of  atoms  in  the  atmosphere  to  those  on  the 
spacecraft  surface  (Bruinsma  and  Biancale,  2003). 

An  increase  in  atmospheric  drag  decreases  orbit  altitude  and  increases  orbit 
velocity.  Thus,  an  inaccurate  drag  forecast  will  result  in  inaccurate  position  forecasts  for 
satellites  in  low-earth  orbit.  In  addition,  increased  drag  over  longer  periods  of  time  will 
decrease  a  satellite’s  operational  lifetime  by  decreasing  its  orbit  altitude  until  it 
experiences  re-entry  (Owens,  et  al.,  2000).  Accurate  characterization  of  thermospheric 
density  is  necessary  in  the  short  term  for  accurate  position  modeling  and  in  the  long  term 
for  accurate  satellite  lifetime  projections. 
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Thermospheric  Density  Measurements 

Accurate  and  continuous  measurements  of  thermospheric  density  have  become 
readily  available  over  the  last  10  years  from  the  Challenging  Minisatellite  Payload 
(CHAMP)  (Bruinsma,  et  al.,  2004)  and  Gravity  Recovery  and  Climate  Experiment 
(GRACE)  (Tapley,  et  al.,  2004)  satellites.  Densities  are  derived  from  on-board 
accelerometers  that  measure  the  electrostatic  force  needed  to  maintain  a  proof  mass  at  the 
center  of  a  cage  located  within  2  mm  of  the  spacecraft’s  center  of  mass.  Since  the 
spacecraft  and  the  proof  mass  respond  to  gravity  in  the  same  way,  the  changes  in  the 
electrostatic  force  needed  to  maintain  the  proof  mass’s  position  reflect  the  spacecraft’s 
response  to  non-gravitational  forces  such  as  thermospheric  drag  (Bruinsma  and  Biancale, 
2003).  The  availability  of  reliable  in-situ  thermospheric  density  measurements  allows 
relevant  comparisons  with  current  modeled  densities  as  well  as  “ground  truth”  data  with 
which  to  test  new  methods  of  modeling  the  thermospheric  environment. 

Thermosphere  as  a  Driven-Dissipative  Thermodynamic  System 

One  approach  for  modeling  the  thermosphere  was  developed  by  Burke  et  al.,  2009 

in  which  the  thermosphere  is  assumed  to  be  a  driven-dissipative  thermodynamic  system. 

The  term  “driven-dissipative”  simply  describes  the  behavior  of  a  system  which  gains 

energy  from  an  input  source,  or  “driver”,  but  then  contains  a  mechanism  which  dissipates 

the  excess  energy  once  the  driver  is  lessened.  This  type  of  system  is  described  by  a 

differential  equation  of  the  same  form  as  that  governing  the  behavior  of  the  disturbance 

storm  time  index  (Dst),  an  index  that  monitors  geomagnetic  activity  at  low  latitudes.  The 

driven-dissipative  approach  uses  empirical  coupling  and  relaxation  constants  to  model  the 
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input  of  energy  to  the  thermosphere  from  the  solar  wind  during  geomagnetic  storms  and 
the  recovery  of  the  thermosphere  back  to  quiet  conditions  after  the  storm  period, 
respectively.  Values  for  these  empirical  constants  were  determined  by  Burke  et  al.,  2009 
by  examining  just  two  storm  periods  during  2004.  Similar  differential  equations  and 
constants  can  be  used  to  model  thermospheric  energy,  exospheric  temperature,  and  Dst. 
Many  existing  thermospheric  density  models  use  exospheric  temperature  as  a  key  input. 
By  obtaining  a  predicted  value  of  exospheric  temperature  from  solar  wind  data,  this 
approach  seeks  to  provide  a  more  accurate  input  for  existing  density  models  that  can  be 
linked  to  solar  wind  models  to  provide  improved  forecast  capabilities. 

Problem 

While  the  driven-dissipative  model  approach  of  Burke  et  al.,  2009  showed 
promising  results  when  compared  to  observed  data  from  GRACE,  it  was  not  applied  to  a 
large  enough  sample  of  storm  events  to  establish  its  general  applicability.  In  later  work 
Burke,  2011  used  the  driven-dissipative  model  to  establish  coupling  constants  for  38 
geomagnetic  storms  between  2002  and  2008.  Burke’s  approach  leaves  several  areas  open 
to  improvement.  This  thesis  expands  on  the  approach  of  Burke  et  al.,  2009  in  the 
following  main  areas: 

1 .  Burke  used  two  storms  in  2004  to  determine  the  value  of  the  relaxation 
constant  and  did  not  allow  it  to  vary  for  other  storms.  This  value  is  suspect  because  Burke 
et.  al,  2009  used  an  early  version  of  GRACE  data  that  has  been  replaced  by  a  revised 
calibration  (Burke,  2011)  (Sutton,  2011).  It  is  also  likely  that  different  storms  will  have 
different  optimal  relaxation  constants.  In  addition,  Burke  used  “trial  and  error 
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comparisons”  (Burke,  et  al.,  2009)  to  determine  values  for  the  coupling  constant  by 
attempting  to  generally  align  model  results  with  the  peak  values  in  observed  data.  A  more 
rigorous  approach  to  determine  the  optimal  values  for  both  the  coupling  and  relaxation 
constants  is  applied  here. 

2.  Burke  et  al.,  2009  treated  energy  input  to  the  thermosphere  from  solar  extreme 
ultra-violet  (EUV)  irradiance  as  a  constant  through  each  storm  period.  In  this  thesis,  it  is 
allowed  to  vary. 

3.  Burke  used  a  simplified  method  of  calculating  observed  orbit-averaged 
GRACE  densities  and  exospheric  temperatures  (Burke,  et  al.,  2009).  In  Burke’s  approach 
orbit-averaged  values  of  density  and  height  were  calculated  from  raw  GRACE 
measurements  and  then  a  quadratic  fit  to  the  Jacchia,  1977  model  (J77)  was  applied  to 
determine  an  orbit-averaged  exospheric  temperature.  This  thesis  modifies  the  orbit¬ 
averaging  technique  and  uses  a  different  application  of  J77  to  produce  more  physically 
accurate  temperatures. 

4.  Burke  modeled  exospheric  temperature  but  current  thermospheric  models  use  a 
global  temperature  parameter  to  model  the  EUV  contribution  to  the  thermospheric  energy 
budget.  J77  uses  a  parameter  known  as  the  arithmetic  mean  temperature,  T\n.  This  thesis 
modifies  Burke’s  approach  to  model  the  arithmetic  mean  temperature. 

Overview 

By  modifying  Burke’s  approach,  this  thesis  provides  a  more  rigorous  test  of  the 
applicability  of  the  driven-dissipative  system  model.  The  result  is  a  more  accurate, 
generalized  model  of  thermospheric  temperatures  using  solar  wind  inputs  as  a  driver. 
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Since  exospheric  temperature  is  used  as  a  parameter  in  existing  thermospheric  models  to 
determine  densities  (Wise,  et  al.,  2012),  a  more  accurate  specification  of  exospheric 
temperature  can  be  used  to  improve  density  forecasts. 

The  remainder  of  this  thesis  is  structured  as  follows.  Section  II  provides 
background  information  on  the  thermospheric  energy  budget,  thermospheric  variability, 
current  thermospheric  models,  and  Burke’s  driven-dissipative  system  model.  Section  III 
details  the  methodology  used  to  develop  the  model  formulation  of  this  thesis.  Section  IV 
presents  the  results  of  the  updated  model  formulation  and  where  appropriate  compares 
the  results  with  Burke’s  earlier  work.  Finally  section  V  presents  conclusions  and 
recommendations  for  future  research. 
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II.  Background 


The  Thermosphere 

The  thermosphere  is  generally  defined  as  the  neutral  part  of  the  Earth’s  upper 
atmosphere  from  95  to  1000  km  above  sea  level.  It  is  characterized  by  a  temperature 
profile  that  increases  with  height  in  its  lower  levels  to  a  maximum  constant  value  which 
is  maintained  to  the  top  of  the  atmosphere  (Schunk  and  Nagy,  2009).  The  top  of  the 
thermosphere  is  defined  as  the  altitude  at  which  neutral  densities  become  low  enough  that 
collisions  become  negligible,  the  atmosphere  can  no  longer  be  treated  as  a  fluid  and 
individual  atoms  and  molecules  have  a  realistic  probability  of  escaping  the  atmosphere  all 
together  (Schunk  and  Nagy,  2009).  This  level  is  referred  to  as  the  exobase  and  the 
temperature  at  this  level,  the  exospheric  temperature,  is  a  major  input  for  many  current 
thermospheric  models. 

Thermospheric  Energy  Input 

There  are  three  main  sources  of  energy  input  to  the  thermosphere:  Extreme 
ultraviolet  (EUV,  X  <  175nm)  irradiance  from  the  sun,  joule  heating,  and  particle 
precipitation  (Knipp,  et  al.,  2004).  Figure  1  shows  the  contribution  of  each  energy  input 
over  the  period  of  solar  cycles  21-23  from  1975  through  2003.  The  lower  gray  curve  in 
Figure  1  shows  the  power  input  to  the  thermosphere  from  particle  precipitation,  the  blue 
curve  represents  the  joule  power  input  and  the  upper  red  curve  depicts  power  input  from 
EUV  irradiance.  EUV  irradiance  in  general  dominates  the  day  side  of  the  thermosphere 
and  is  closely  associated  with  the  1 1-year  solar  cycle.  Joule  heating  and  particle 
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precipitation  are  most  important  in  the  auroral  zones  and  are  closely  associated  with 
geomagnetic  activity  (Knipp,  et  al.,  2004).  Each  of  the  three  energy  sources  is  discussed 
in  the  subsequent  sections. 
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Figure  1:  Power  input  to  the  thermosphere  by  particle  precipitation  (bottom  gray  line),  joule  heating 
(blue  line),  and  EUV  irradiance  (red  line)  for  each  day  from  1975  to  2003.  Adapted  from  Knipp,  et  al., 

2004. 


Solar  EUV  Irradiance 

Solar  EUV  irradiance,  emitted  from  the  sun’s  chromosphere  and  corona,  is 
usually  the  dominant  contributor  to  thermospheric  energy.  From  1975-2003  solar 
irradiance  made  up  an  average  of  78%  of  the  total  energy  input  to  the  thermosphere 
(Knipp,  et  al.,  2004).  The  energy  is  deposited  mainly  in  the  layer  from  150-200  km 
(Knipp,  et  al.,  2004)  via  absorption  by  neutrals,  primarily  O2  and  N2  (Schunk  and  Nagy, 
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2009).  Because  a  portion  of  the  absorbed  energy  goes  into  dissociation  and  ionization,  the 
heating  efficiency  is  limited  to  around  50%  (Knipp,  et  al.,  2004).  The  solar  irradiance 
contribution  to  thermospheric  energy  varies  by  100%  or  more  over  the  course  of  a  solar 
cycle  as  seen  in  Figure  1.  The  day  to  day  variation  is  much  smaller  during  solar 
minimum  than  near  solar  maximum. 

Joule  Heating 

Joule  heating  is  the  process  in  which  an  electric  current  passes  through  the 
thermosphere  resulting  in  resistance  and  heating  of  the  neutral  gas  (Qian  and  Solomon, 
2011).  On  average,  joule  heating  accounts  for  16%  of  the  total  energy  input  to  the 
thermosphere  via  deposition  mainly  from  1 10-140  km  (Knipp,  et  al.,  2004).  The  energy 
source  for  thermospheric  joule  heating  is  the  solar  wind  which  interacts  with  the 
magnetosphere  to  create  electric  fields  that  map  into  the  thermosphere  and  drive  currents. 
Lu  et  al.,  1998  showed  that  on  average  about  60%  of  the  solar  wind  energy  that  is 
transferred  to  the  magnetosphere  is  deposited  in  the  thermosphere.  During  storm  times 
the  amount  of  solar  wind  energy  deposited  in  the  thermosphere  through  the 
magnetosphere  can  reach  80%  (Lu,  et  al.,  1998).  Since  joule  power  input  is  caused  by 
currents,  it  can  be  monitored  with  indices  that  respond  to  ionospheric  or  magnetospheric 
currents  such  as  the  AE  index,  which  monitors  the  auroral  electrojet,  and  the  Dst  index, 
which  monitors  the  ring  current  (Knipp,  et  al.,  2004). 

While  joule  power  input  is  generally  much  smaller  than  solar  irradiance,  it 
exhibits  more  variability.  Since  joule  heating  is  over  90%  efficient  in  transferring  power 
to  the  thermosphere  (compared  to  50%  efficiency  for  the  solar  and  particle  inputs)  any 
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change  in  available  power  is  readily  transferred  to  the  thermosphere  (Knipp,  et  al.,  2004; 
Thayer  and  Semeter,  2004).  When  looking  at  the  100  days  with  the  highest  total  power 
inputs  of  from  solar  cycles  21  -  23,  Knipp,  et  al.,  2004  found  that  solar  irradiance 
increased  50%  above  its  average  value  while  joule  power  increased  by  over  600%  above 
its  average.  During  large  geomagnetic  storms  the  joule  power  input  becomes  the 
dominant  power  source  for  the  thermosphere  and  when  combined  with  the  particle 
precipitation  power  input,  the  power  input  due  to  geomagnetic  activity  accounts  for  65% 
of  the  total  (Knipp,  et  al.,  2004).  When  examining  thermospheric  variability  on  short  time 
scales  joule  power  becomes  the  most  important  term. 

Particle  Precipitation 

Another  way  that  energy  is  transferred  from  the  solar  wind  to  the  thermosphere  is 
via  precipitation  of  electrons.  Solar  wind  electrons  travel  along  open  magnetic  field  lines 
or  through  the  magnetotail  into  the  auroral  zone  where  they  are  absorbed  (Prolss,  2004), 
primarily  from  100-120  km  (Knipp,  et  al.,  2004).  Since  some  of  the  electron  energy  goes 
into  ionization,  rotational,  or  vibrational  states  the  heating  efficiency  for  particle 
precipitation  is  limited  to  around  50%.  On  average,  particle  precipitation  accounts  for  6% 
of  the  total  power  input  to  the  thermosphere  (Knipp,  et  al.,  2004).  During  geomagnetic 
activity,  the  magnetosphere  interacts  with  the  solar  wind  magnetic  field  resulting  in  more 
open  field  lines  and  more  available  paths  for  electrons  to  reach  and  transfer  power  to  the 
thermosphere  (Prolss,  2004).  Strong  geomagnetic  storms  result  in  an  increase  in  power 
input  due  to  particle  precipitation  of  up  to  200%  compared  with  average  values  (Knipp,  et 
al.,  2004). 
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Thermospheric  Energy  Loss 

On  long  time  scales  the  thermospheric  system  is  in  a  state  of  equilibrium  where 
the  energy  input  is  equal  to  the  energy  loss  as  evidenced  by  observations  that  show 
thermospheric  temperatures  do  not  increase  or  decrease  indefinitely.  One  of  the  major 
loss  mechanisms  for  thermospheric  energy  is  emission  by  nitric  oxide  (NO)  at  5.3  pm 
(Sharma,  et  ah,  1996).  Radiation  at  5.3  pm  is  not  readily  absorbed  by  any  major 
atmospheric  constituent  so  energy  at  this  wavelength  is  able  to  escape  into  space. 

To  maintain  equilibrium,  there  must  be  a  mechanism  during  storm  time  by  which 
the  excess  energy  input  to  the  thermosphere  via  joule  heating  and  particle  precipitation  is 
dissipated  as  the  thermosphere  relaxes  to  its  pre-storm  state.  Since  the  production  rate  of 
NO  is  highly  dependent  on  temperature  (Bailey,  et  al.,  2002),  the  high  thermospheric 
temperature  during  geomagnetic  storms  leads  to  increased  NO  densities  resulting  in 
increased  cooling  rates.  Mlynczak  et  al.,  2005  found  that  increased  NO  emissions  during 
geomagnetic  storming  accounted  for  roughly  94%  of  the  added  thermospheric  energy 
loss  during  the  recovery  period.  The  rest  of  the  energy  loss  increase  can  be  accounted  for 
by  increased  CO2  emissions  at  1 5  pm  (2%)  and  increased  conduction  between  the 
thermosphere  and  mesosphere  (4%)  (Mlynczak,  et  al.,  2005). 

Thermospheric  Variability 

Solar  EUV  irradiance  is  the  primary  energy  input  to  the  thermosphere,  while  joule 

heating  and  particle  precipitation  are  secondary  the  majority  of  the  time.  Both  the  primary 

and  secondary  drivers  result  in  variability  in  thermospheric  densities  and  temperatures  on 

different  time  scales  and  each  can  be  accounted  for  through  the  use  of  various 
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observations  and  proxies.  Variations  on  specific  temporal  and  spatial  time  scales  will  be 
discussed  below. 

Solar  Cycle  Variability 

The  sun  exhibits  a  cycle  between  solar  minimum  and  solar  maximum  with  a 
period  of  roughly  1 1  years,  characterized  in  part  by  changes  in  solar  irradiance  (Figure  1). 
During  solar  maximum  there  are  many  more  active  regions  on  the  sun  resulting  in 
increased  irradiance,  increased  flaring,  and  more  frequent  coronal  mass  ejections  which 
in  turn  increase  geomagnetic  activity.  This  periodic  irradiance  variation,  along  with  the 
increase  in  geomagnetic  activity  as  a  lesser  factor,  generates  a  similar  variation  in 
exospheric  temperature,  thermospheric  energy  and  density  at  the  earth.  The 
thermospheric  density  at  a  given  altitude  during  solar  max  can  be  up  to  ten  times  more 
than  the  density  at  that  same  altitude  during  solar  min  (Qian  and  Solomon,  2011). 

Semiannual  Variability 

Thermospheric  density  varies  on  a  semiannual  basis  with  maximums  at  the 
equinoxes  and  minimums  near  the  solstices.  This  variation  was  first  identified  by 
Paetzold  and  Zchomer,  1961  when  they  showed  that  the  difference  between  minimum 
and  maximum  is  more  than  100%.  Semi-annual  variability  is  driven  primarily  by  the 
variation  in  the  distance  from  the  sun  to  the  earth  which  causes  differences  in  solar 
irradiation.  Between  this  variation  and  the  solar  cycle  variation  described  above,  it  is 
clear  that  even  with  geomagnetic  activity  removed  from  consideration  the  density  of  the 
thermosphere  fluctuates.  Any  attempt  to  model  densities  accurately  must  account  for 
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variations  in  both  the  solar  and  geomagnetic  contributions  if  it  is  to  accurately 
characterize  the  thermospheric  environment. 

Solar-Rotation  Variability 

The  sun  rotates  differentially  with  an  average  period  of  27  days  and  during  this 
rotational  period  active  regions  of  the  sun  appear  and  disappear  from  the  Earth’s  view. 
Since  active  regions  can  persist  for  several  months  they  may  come  into  and  disappear 
from  the  Earth’s  view  multiple  times  during  their  lifetime.  Active  regions  are  associated 
with  increased  solar  irradiance  and  geomagnetic  activity  and  therefore  solar  rotation 
results  in  periodic  changes  in  irradiance  and  geomagnetic  activity.  This  periodic 
variability  in  irradiance  and  geomagnetic  activity  results  in  a  variability  of  up  to  100%  in 
thermospheric  density  during  solar  maximum  (Qian  and  Solomon,  2011). 

Multi-Day  Variability 

Variations  in  the  solar  wind  caused  by  high-speed  streams  (HSS)  result  in  low 
levels  of  geomagnetic  activity  and  can  therefore  impact  the  thermosphere  via  increased 
joule  heating  and  particle  precipitation.  Observations  during  the  declining  phase  of  solar 
cycle  23  showed  periodic  variations  in  the  source  of  HSS,  coronal  holes  (Temmer,  et  al., 
2007).  Similarly  Lei  et  al.,  2008,  found  a  9-day  periodic  variation  in  neutral  density 
observations  from  the  CHAMP  satellite  in  2005.  The  magnitude  of  these  variations  is 
smaller  than  those  due  to  solar  rotation,  roughly  30  -  50%  in  density. 

Diurnal  Variability 

As  expected,  the  large  disparity  in  solar  irradiance  between  the  day  and  night 
sides  of  the  thermosphere  results  in  a  large  density  variation  between  the  two.  Mueller  et 
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al.,  2009,  found  that  during  geomagnetic  quiet  periods  the  density  on  the  day  side  was 
roughly  twice  that  on  the  night  side.  By  taking  orbit  averages  of  density  measurements 
from  polar  orbiting  satellites  such  as  GRACE  the  diurnal  variation  can  be  effectively 
averaged  out  of  observed  data. 

Short  Term  Variability 

Density  variations  on  time  scales  of  minutes  to  hours  can  be  caused  by  rapidly 
changing  energy  inputs  to  the  system  from  solar  flares  or  geomagnetic  storms  associated 
with  coronal  mass  ejections  (CMEs)  or  high-speed  streams.  Solar  flares  cause  rapid 
increases  in  EUV  and  X-ray  irradiance  leading  to  heating  and  expansion  of  the  upper 
atmosphere  (Pawlowski  and  Ridley,  2008).  Thermospheric  density  increases  depend  on 
the  flare’s  intensity,  location,  and  the  details  of  the  flare’s  spectral  enhancement.  Density 
increases  of  up  to  40%  have  been  observed  in  response  to  long  duration  ( >  40  min)  X- 
class  flares  (Qian  and  Solomon,  2011). 

Geomagnetic  storms  also  result  in  increased  energy  inputs  to  the  thermosphere 
however  the  process  by  which  the  energy  is  deposited  is  different.  During  geomagnetic 
storms  energy  is  transferred  from  the  solar  wind  to  the  thermosphere  via  joule  heating 
and  particle  precipitation  at  high  (auroral)  latitudes.  Joule  heating  in  the  thermosphere  is 
the  dominant  form  of  energy  transfer  over  particle  precipitation  during  geomagnetic 
storms  (Wilson,  et  al.,  2006).  The  energy  deposited  at  high  latitudes  is  propagated 
throughout  the  thermosphere  via  circulation  and  atmospheric  gravity  waves  over  a  time 
period  of  several  hours  (Bruimsma,  et  al.,  2006).  The  focus  of  this  research  is  to  better 
characterize  the  impact  of  geomagnetic  storms  on  the  thermosphere. 
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Geomagnetic  Storming 


A  geomagnetic  storm  has  been  defined  by  Prolss,  2004  as  “an  event  of  strongly 
enhanced  dissipation  of  solar  wind  energy  in  the  near-Earth  space  environment.”  During 
geomagnetic  storming  both  the  joule  heating  and  particle  precipitation  energy  inputs  to 
the  thermosphere  are  enhanced.  The  dominant  factors  in  determining  the  amount  of 
energy  transfer,  and  therefore  the  strength  of  a  geomagnetic  storming  event,  is  the 
component  of  the  interplanetary  magnetic  field  in  the  z  direction,  Bz  and  the  length  of 
time  Bz  is  in  the  negative  z  direction.  The  z  direction  is  defined  by  Geocentric  Solar 
Magnetospheric  (GSM)  coordinates  shown  in  Figure  2. 


Figure  2:  Illustration  of  the  geocentric  solar  magnetospheric  (GSM)  coordinate 
system.  The  origin  is  the  center  of  the  earth,  the  x  axis  points  toward  the  sun,  the 
y  axis  is  perpendicular  to  both  x  and  the  geomagnetic  dipole  axis,  and  the  z  axis 
completes  the  set  with  positive  pointing  north.  (Knecht  and  Shuman,  1985) 


Solar  wind  energy  is  transferred  to  the  thermosphere  through  the  magnetosphere 
via  a  dynamo  of  conductive  solar  wind  plasma  moving  across  the  Earth’s  magnetic  field 
lines.  This  dynamo  is  made  possible  by  an  “open  magnetosphere”  magnetic  field 
configuration  created  by  the  interaction  between  a  southward  interplanetary  magnetic 
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field  and  the  Earth’s  dipole  magnetic  field  shown  in  Figure  3.  The  southward  Bz  interacts 
with  Earth’s  northward-pointing  magnetic  field,  weakening  the  field  on  the  day  side  of 
Earth  and  resulting  in  an  increased  number  of  open  magnetic  field  lines  (Prolss,  2004). 

An  open  magnetic  field  line  has  one  footpoint  on  Earth  in  the  auroral  region  and  the  other 
in  space  (Prolss,  2004).  These  open  magnetic  field  lines  provide  pathways  that  allow 
energetic  particles  to  reach  the  thermosphere,  increasing  power  input  from  particle 
precipitation,  and  allow  an  electric  dynamo  to  transfer  energy  from  the  solar  wind  to  the 
thermosphere  via  joule  heating. 
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Figure  3:  The  interaction  between  a  southward  Interplanetary  Magnetic  Field  (Bz  south)  and  the 
Earth’s  dipole  magnetic  field  is  shown.  The  result  is  open  magnetic  field  lines,  with  one  footpoint  near 
the  polar  cap  and  the  other  in  interplanetary  space.  This  configuration  is  referred  to  as  the  open 
magnetosphere.  Figure  from  Prolss,  2004. 


Prolss,  2004  describes  the  energy  transfer  process  as  follows.  With  an  open 
magnetosphere  the  Earth’s  magnetic  field  lines  originating  near  the  polar  cap  are  not 
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closed  but  open  extending  into  the  interplanetary  medium,  as  shown  in  Figure  4.  As  the 
solar  wind  flows  across  this  magnetic  field  the  charged  particles  experience  a  Lorentz 


Interplanetary  space 


Figure  4:  The  interaction  between  the  solar  wind  and  the  open  magnetosphere 
configuration  is  shown.  Figure  from  Prolss,  2004. 


force,  causing  the  positively  charged  particles  to  be  deflected  towards  the  dawn  side  and 
the  negative  particles  to  be  deflected  towards  the  dusk  side.  The  resultant  charge 
separation  creates  a  polarization  electric  field,  eP,  which  builds  up  until  the  force  on 
charged  particles  due  to  the  polarization  field  matches  that  due  to  the  Lorentz  force,  as 
shown  in  Equation  (2) 
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nsqs£p+nsqsuswY.Bz=0 


(2) 


where  ns  is  the  number  density  of  a  given  species,  qs  is  the  charge  of  a  given  species,  and 
usw  is  the  solar  wind  velocity.  Solving,  we  see  that  the  polarization  electric  field  is  equal 
to  the  negative  solar  wind  velocity  crossed  with  the  z  component  of  the  magnetic  field. 
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This  quantity  is  also  known  as  the  electric  dynamo  field,  8dyn,  and  is  shown  in  Figure  4. 
The  dawn  to  dusk  electric  dynamo  field  maps  along  the  magnetic  field  lines  to  the  polar 
cap  region  where  it  drives  a  current,  denoted  in  Figure  4  as  jp,  that  deposits  energy  into 
the  thermosphere/ionosphere  system  via  joule  heating. 

The  energy  input  to  the  thermosphere  by  joule  heating  is  extracted  from  the  solar 
wind  and  manifested  through  a  reduction  in  solar  wind  velocity.  The  electric  dynamo 
field  drives  a  current  in  the  magnetosphere,  denoted  in  Figure  4  by  jdyn,  which  interacts 
with  the  interplanetary  magnetic  field  to  produce  a  force  in  the  direction  opposing  the 
solar  wind  flow  and  decreasing  the  flow  velocity. 


Fr  =  1,  xB  (4) 

B  J  dyn  z 

The  current  loop  between  the  polar  cap  current,  jp,  and  the  dynamo  current,  jdyn,  is  closed 
by  the  region  one  Birkeland  currents,  js,  shown  in  Figure  4.  Region  one  currents  are 
defined  as  currents  originating  on  the  poleward  boundary  of  the  auroral  oval  (Prolss, 
2004). 
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An  open  magnetosphere  configuration  is  necessary  to  create  the  enhanced  joule 
heating  and  particle  precipitation  power  input  to  the  thermosphere  observed  during 
geomagnetic  storming.  The  two  main  solar  phenomena  which  lead  to  strong  southward 
Bz,  creating  an  open  magnetosphere  and  geomagnetic  storming,  are  coronal  mass 
ejections  (CMEs)  and  co-rotating  interaction  regions  (CIRs). 

Coronal  Mass  Ejections 

A  CME  is  a  large  emission  of  mass  from  the  sun,  on  the  order  of  1012  -  1013  kg,  at 
speeds  of  50  -  1800  km/s  with  an  average  kinetic  energy  ranging  from  10  to  10  J 
(Prolss,  2004).  CMEs  are  accelerated  outward  from  the  sun  by  magnetic  forces  in  the 
sun’s  corona.  Depending  on  the  orientation  of  the  magnetic  field  within  the  ejected 
material,  a  CME’s  encounter  with  earth  can  produce  a  southward  Bz  along  with  enhanced 
solar  wind  velocity  and  density  resulting  in  geomagnetic  storming  (Prolss,  2004). 

Co-Rotating  Interaction  Regions 

CIRs  have  their  source  on  the  sun  at  the  boundaries  between  coronal  holes  and 
coronal  streamers.  Coronal  holes  are  a  source  of  high-speed  solar  wind  streams  while 
coronal  streamers  are  a  source  of  low-speed  solar  wind  flow  (Prolss,  2004).  As  the  solar 
wind  propagates  out  from  the  sun  towards  Earth  the  difference  in  velocity  between  the 
two  regions  results  in  a  compression  of  the  solar  wind  plasma  in  the  area  where  the  high¬ 
speed  stream  interacts  with  the  low-speed  flow.  This  area  of  compression  is  defined  as  a 
CIR.  When  solar  plasma  leaves  the  sun  as  the  solar  wind  it  carries  with  it  a  “frozen-in” 
magnetic  field  with  the  same  orientation  as  its  source  region  on  the  sun.  The  magnetic 
field  is  compressed  along  with  the  plasma  inside  the  CIR.  If  the  frozen-in  magnetic  field 
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was  already  oriented  southwards,  the  amplification  inside  the  CIR  is  sufficient  to  produce 
geomagnetic  storming  when  the  CIR  encounters  Earth’s  magnetic  field  (Prdlss,  2004). 

Storm  Type  Characteristics 

A  geomagnetic  storm  produced  by  a  CME  is  distinct  from  one  produced  by  a  CIR 
in  several  ways  (Borovsky  and  Denton,  2006).  The  rate  of  CME-driven  storm  occurrence 
peaks  during  solar  maximum  and  is  smallest  during  solar  minimum  (Webb,  1991)  while 
the  frequency  of  CIR-driven  storms  is  the  highest  during  the  declining  phase  of  the  solar 
cycle  (Mursula  and  Zeiger,  1996).  The  occurrence  pattern  for  CME-driven  storms  is 
irregular  with  no  characteristic  spacing  between  events  while  CIR-driven  storms  are 
characterized  by  a  27-day  periodicity  due  to  the  rotation  of  their  source  regions,  coronal 
holes,  on  the  sun  (Borovsky  and  Denton,  2006).  CME-driven  storms  are  more  effective 
than  CIR-driven  storms  in  producing  highly  negative  Dst  values  (Dst  <  -100  nT)  and  are 
usually  characterized  by  a  shock  in  the  solar  wind  flow,  evidenced  by  a  sharp  increase  in 
solar  wind  velocity  and  density  (Borovsky  and  Denton,  2006).  CIR-driven  storms 
normally  produce  less  extreme  Dst  values  and  exhibit  a  more  gradual  commencement. 
These  differences  were  used  to  classify  the  storms  used  in  this  thesis  as  either  CME  or 
CIR  storms. 

Figure  5  shows  typical  solar  wind  profiles  for  both  CME  (top)  and  CIR  (bottom) 
storms.  The  CME  storm  has  an  extreme  Dst  minima  of  -181  nT  while  the  CIR  storm  does 
not  drop  below  -50  nT.  The  start  of  the  CME  storm  is  evident  in  the  rapid  rise  in  solar 
wind  pressure  and  velocity  around  JD  250.7.  In  contrast,  the  CIR  storm  exhibits  a  gradual 
increase  in  solar  wind  pressure  and  veolcity  between  JD  191.5  and  192. 
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Figure  5:  Typical  solar  wind  signatures  resulting  from  a  coronal  mass  ejection 
(CME)  driven  geomagnetic  storm  (top)  and  a  co-rotating  interaction  region  (CIR) 
driven  geomagnetic  storm  (bottom).  From  top  to  bottom  the  plots  show  Dst  in 
nano-Tesla,  the  z-component  of  the  interplanetary  magnetic  field  in  nano-Tesla 
(GSM  coordinates),  solar  wind  pressure  (P)  in  nano-pascals,  solar  wind  velocity 
(V)  in  km/s,  and  the  resulting  magnetospheric  electric  field  value,  in 
milivolts/meter  as  functions  of  julian  date  (JD)  counted  as  days  since  1  January  of 
the  given  year. 
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Thermospheric  Driver  Proxies 

Since  both  the  primary  (solar  EUV  irradiance)  and  secondary  (geomagnetic 
activity)  sources  of  thermospheric  energy  have  historically  been  difficult  to  measure 
directly  various  proxies  and  indices  are  used  to  quantify  their  variation  for  use  in  models. 
In  some  thermospheric  models  geomagnetic  activity  has  been  accounted  for  through  the 
use  of  the  ap  index.  The  ap  index  is  a  linear  index  ranging  from  zero  to  400  that  is  derived 
using  the  deviation  from  the  standard  magnetic  field  values  measured  at  13  locations 
worldwide  at  geomagnetic  latitudes  ranging  from  42  to  62  degrees  (Helmholtz  Centre 
Potsdam  GFZ,  2012).  Values  are  computed  every  three  hours  for  ap  and  daily  averages 
are  computed  and  reported  as  Ap.  The  fact  that  the  ap  index  is  measured  at  mid  latitudes 
results  in  a  failure  to  detect  the  full  impact  of  large  geomagnetic  storms  (Bowman,  et  al., 
2008)  due  to  distortion  from  the  equatorward  movement  of  the  auroral  electrojet. 

Another  measure  of  geomagnetic  activity  is  the  disturbance  storm  time  index 
(Dst).  Dst  is  measured  hourly  at  four  different  near-equatorial  observatories  and  it 
measures  the  variations  in  the  Earth’s  magnetic  field  resulting  from  changes  in  the 
magnetospheric  ring  current.  Since  the  ring  current  responds  directly  to  energy  inputs 
from  the  solar  wind,  it  is  enhanced  during  periods  of  geomagnetic  storming.  Dst  is 
measured  in  nano-Tesla  (nT)  and  during  quiet  conditions  it  is  usually  near  zero.  Storming 
conditions  are  indicated  by  negative  values  and  the  more  negative  the  value  the  stronger 
the  storm.  Because  of  the  equatorial  location  of  its  observation  stations,  Dst  is  not 
influenced  by  the  auroral  zone  and  is  able  to  detect  more  fully  the  energy  enhancements 
to  the  ring  current  caused  by  strong  geomagnetic  storms.  The  Dst  index  has  been  adopted 
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for  use  in  some  recent  thermospheric  models  such  as  Jacchia-Bowman  2008  (Bowman,  et 
al„  2008). 

TheF  10.7  index  has  long  been  the  standard  proxy  for  EUV  flux.  Since  the 
atmosphere  absorbs  virtually  all  of  the  EUY  radiation  before  reaches  the  surface  it  is  not 
possible  to  measure  EUV  flux  at  a  surface  based  observatory.  Instead,  EUV  flux  values 
can  be  inferred  using  measurements  of  the  solar  radio  flux  at  a  wavelength  of  10.7  cm  at 
the  Earth’s  surface.  This  10.7  cm  flux  has  been  shown  to  correlate  well  with  actual  EUV 
flux.  F  io.7  values  are  observed  at  the  Pentictin  Radio  Observatory  in  British  Columbia, 
Canada  daily  at  2000Z  (local  noon).  Daily  F10.7  values,  along  with  a  longer  term  81  or 
162-day  average,  have  been  used  in  many  models  to  account  for  the  variation  in  EUV 
flux  (Tascione,  1994).  Unfortunately,  the  observed  nature  of  the  F10.7  index  and  its  once- 
daily  time  resolution  have  limited  models  making  use  of  it  as  a  input. 

Partly  in  an  effort  to  overcome  these  limitations,  the  first  full-spectrum  solar 
irradiance  model,  SOLAR2000,  was  developed  by  Tobiska  et  al.  in  2000.  This  model 
includes  a  new  EUV  proxy  index,  E10.7,  which  is  in  the  same  units  of  the  standard  F10.7 
index  so  as  to  enable  its  use  in  existing  modeling  applications.  The  E10.7  has  several 
advantages  over  the  F10.7  including  the  availability  of  high  temporal  resolution  data  rather 
than  the  once-daily  F10.7  and  the  ability  to  forecast  E10.7  values  into  the  future  which  does 
not  exist  with  the  observed  F10.7  index.  Some  recent  models,  such  as  HASDM,  have 
adopted  E10.7  to  replace  F10.7  for  these  and  other  reasons  (Storz,  et  al.,  2005). 

The  E10.7  models  total  integrated  EUV  emissions  from  both  the  chromosphere  and 
the  corona  while  the  F10.7  proxy  only  captures  coronal  emissions.  By  providing  a  more 
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complete  picture  of  total  EUV  irradiance  the  E10.7  is  a  more  representative  proxy  for  the 
impact  of  EUV  irradiance  on  the  thermosphere.  However,  it  leads  to  differences  when 
compared  with  the  longtime-standard  F10.7.  Tobiska  et  al.,  2000  found  that  F10.7  exhibited 
more  variability  than  E10.7,  as  much  as  +/-  20%  during  comparisons  ran  for  July,  1982. 
The  increased  variability  of  the  F10.7  was  due  to  the  fact  that  it  does  not  measure 
chromospheric  emissions,  which  tend  to  smooth  out  the  E10.7  values. 

Thermospheric  Models 

These  indices  and  proxies,  along  with  historic  and  real-time  observations,  have 
been  used  to  create  many  different  models  of  the  thermosphere.  The  following  sections 
briefly  describe  three  thermospheric  models  relevant  to  this  thesis. 

Jacchia  Models 

Jacchia  developed  a  model  of  the  thermosphere  in  1970  (J70)  (Jacchia,  1970)  and 
an  updated  version  in  1977  (J77)  (Jacchia,  1977)  that  are  still  used  as  a  baseline  today. 
The  Jacchia  models  are  static  models  which  were  developed  using  thermospheric 
densities  calculated  from  satellite  drag  and  mass  spectrometer  measurements.  They  are 
based  on  the  assumption  that  the  thermosphere  is  in  thermal  diffusion  equilibrium, 
meaning  that  the  heat  inputs  to  the  thermosphere  equal  heat  losses.  The  J77  model 
assumes  the  mesopause,  the  bottom  of  the  thermosphere,  is  at  an  altitude  of  90  km  with  a 
temperature  of  188K  and  a  mass  density  of  3.43  Model  temperatures  rise  as  a 
function  of  altitude  from  the  minimum  value  at  90km,  pass  through  an  inflection  point  at 
125km,  and  increase  asymptotically  to  the  given  exospheric  temperature,  Too .  Too 
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uniquely  defines  the  temperature  profile.  Once  the  temperature  profile  is  determined, 
densities  are  calculated  by  integrating  the  thermal  diffusion  equation,  Equation  (5), 


dn  m.g  dT ,  \ 

n.  RT  T  V  lJ 


(5) 


where  the  index  i  denotes  the  ith  species,  n  is  the  number  density,  m  is  the  mass,  g  is 
gravity,  a  is  the  thermal  diffusion  coefficient,  T  is  the  temperature  and  R*  is  the  universal 
gas  constant.  The  J77  model  includes  six  species:  N2,  O2,  O,  Ar,  He,  and  H.  The  total 
mass  density  at  a  given  altitude  can  be  calculated  by  simply  summing  the  product  n  jn, 
over  all  species  (Wise,  et  al.,  2012).  Through  this  process,  tables  are  produced  that  give 
density  profiles  for  a  given  exospheric  temperature  input. 

Variations  due  to  solar  changes  and  geomagnetic  activity  are  accounted  for  in  the 
Jacchia  models  either  solely  through  perturbations  to  the  temperature  profile  (J70)  or 
through  perturbations  to  both  the  temperature  and  resulting  density  profiles  (J77).  The 
J77  model  accounts  for  variations  in  EUV  energy  input  by  using  the  F10.7  proxy  and  an 
F10.7  index  value  averaged  over  six  solar  rotations  (162-days),  F  10.7a,  to  compute  a 
geomagnetic-quiet  (defined  as  Ap  =  0)  arithmetic-mean  exospheric  temperature,  T i/2uv- 
The  arithmetic  mean  temperature,  T 1/2,  is  defined  as  the  average  of  the  nighttime 
minimum  exospheric  temperature,  To,  and  the  daytime  maximum  exospheric  temperature, 
Tm,  which  occur  in  opposite  hemispheres  at  0524  and  1648  Local  Standard  Time  (LST), 
respectively  (Jacchia,  1977).  T1/2  is  related  to  Too  at  any  given  location  via  a  conversion 
factor  dependent  on  latitude,  local  time  and  solar  declination  angle.  Using  the  J77  model, 
unique  temperature  and  density  profiles  can  be  computed  for  any  location  given  T1/2.  The 
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J70  model  uses  a  similar  process  but  instead  of  T\a  its  global  temperature  parameter  is 
To.  The  tables  in  the  J70  and  J77  models  form  the  basis  of  many  current  thermospheric 
models.  For  this  thesis,  J77  serves  as  the  link  between  observed  neutral  density 
measurements  from  GRACE  data  and  an  “observed”  exospheric  temperature  used  for 
comparison  with  the  exospheric  temperature  calculated  using  Burke’s  driven-dissipative 
model. 

High  Accuracy  Satellite  Drag  Model  (HASDM) 

The  Jacchia  models  have  been  improved  through  the  years  but  continue  to  be 
limited  by  their  use  of  proxies  to  measure  actual  thermospheric  conditions  as  well  as  their 
reliance  on  a  static  and  limited  set  of  observed  data  upon  which  their  empirical  fits  are 
based.  These  limitations,  along  with  others,  prevent  satellite  position  error  from 
decreasing  below  15%  (Marcos,  et  al.,  2007).  To  address  this  problem  the  Air  Force 
Space  Command  Battlelab  created  HASDM,  the  Air  Force’s  current  operational 
thermospheric  density  model,  in  2004  (Storz,  et  al.,  2005). 

HASDM  makes  use  of  the  ap  index  to  characterize  geomagnetic  activity.  To 
characterize  EUV  flux  HASDM  uses  the  E10.7  index  from  the  SOLAR2000  model 
described  by  Tobiska  et  al.,  2000.  The  critical  advance  of  the  HASDM  approach  is  the 
use  of  near  real-time  observed  density  data.  The  model  uses  data  from  the  observed  drag 
on  a  set  of  about  80  calibration  satellites  to  create  spatially  varying  density  corrections 
every  three  hours.  These  corrections  are  used  in  conjunction  with  a  modified  J70  model 
to  produce  a  global  density  forecast  up  to  72  hours  into  the  future.  This  approach  of 
relying  on  observed  data  in  real  time  to  dynamically  update  and  correct  density 
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predictions  helped  reduce  positional  errors  down  to  5%  for  the  calibration  satellites  and 
down  to  8%  for  all  tracked  objects  during  quiet  conditions  (Storz,  et  ah,  2005). 
Unfortunately,  HASDM  does  not  perform  as  well  during  geomagnetic  storming 
conditions.  During  storm  periods  neutral  density  errors  increase  by  roughly  30%,  from 
13%  during  quiet  conditions  (Ap  =  0)  to  17%  during  storming  conditions  (Ap>  100) 
(Marcos,  et  ah,  2010).  HASDM  leaves  room  for  improved  characterization  of  storming 
conditions. 

Jacchia-Bowman  2008  (JB2008)  Model 

JB2008  is  an  empirical  model  which  uses  density  inputs  from  Air  Force  daily 
density  values  (computed  using  tracking  data  from  around  100  calibration  satellites)  and 
HASDM  as  well  as  CHAMP  and  GRACE  accelerometer  data  (Bowman,  et  ah,  2008). 
JB2008  uses  the  F10.7  index  and  the  81 -day  average  F10.7  index  along  with  26  -  34  nm 
integrated  EUV  flux  data  from  the  Solar  Heliospheric  Observatory  (SOHO)  satellite, 
chromospheric  and  photospheric  active  region  activity  data  measured  by  the  Solar 
Backscatter  Ultraviolet  (SBUV)  spectrometer,  and  X-ray  emission  data  from  GOES  X- 
ray  spectrometers  to  compute  To.  This  approach  allows  the  JB2008  model  to  capture  not 
only  solar  cycle  and  semi-annual  solar  irradiance  variations  but  also  measure  shorter  term 
variations  on  the  scale  of  the  27-day  solar  rotation  period. 

Another  advance  of  the  JB2008  model  is  its  use  of  Dst  to  measure  geomagnetic 
activity  rather  than  the  ap  index  used  by  previous  models.  It  is  a  better  input  to 
thermospheric  models  than  ap  because  ap  responds  mainly  to  ionospheric  currents  rather 
than  magnetospheric  ones.  Since  the  energy  deposited  into  the  thermosphere  during 


27 


geomagnetic  storms  comes  from  the  solar  wind  through  the  magnetosphere  it  is 
reasonable  to  use  an  input  that  primarily  measures  magnetospheric  conditions  like  Dst.  In 
addition,  ap  is  determined  by  observatories  at  latitudes  from  42  to  62  degrees  which  can 
incorrectly  characterize  energy  inputs  during  severe  storms  due  to  the  equatorward 
movement  of  the  auroral  electrojet.  Since  ground-based  observatories  are  immobile, 
significant  electrojet  movement  during  storm  time  leads  to  underestimates  of  storm 
impacts  (Huang  and  Burke,  2004).  Dst  responds  to  the  ring  current  and  is  derived  from 
measurements  at  four  equatorial  observatories  not  impacted  by  auroral  electrojets.  Using 
Dst  as  an  input,  a  geomagnetic  activity  contribution  to  To  is  calculated  and  then  used  to 
generate  a  density  profile. 

Modeling  the  Thermosphere  as  a  Driven-Dissipative  Thermodynamic  System 

While  thermospheric  models  have  made  advances  in  accuracy  recently,  they  are 
still  physically  limited  by  the  lack  of  a  direct  link  between  the  solar  wind  and  the 
thermosphere  which  is  the  dominant  source  of  energy  during  geomagnetic  storming.  The 
driven-dissipative  approach  attempts  to  solve  this  problem  by  linking  the  thermosphere  to 
the  solar  wind  using  the  electric  field  of  the  magnetosphere  as  the  primary  driver  during 
geomagnetic  storm  conditions. 

Burke  et  al.,  2009  used  neutral  density  observations  from  the  GRACE  satellite 

along  with  the  J77  model  (Jacchia,  1977)  to  calculate  thermospheric  energies,  Eth,  as  a 

function  of  time  during  2004.  Magnetospheric  electric  field  magnitudes,  sys,  were 

computed  using  observed  solar  wind  data  and  plotted  as  a  function  of  time  along  with  the 

Eth  data.  Figure  6  (Burke,  et  al.,  2009)  shows  that  Eth  decays  to  pre-disturbance  levels 

28 


when  s vs  drops  to  pre-disturbance  levels  and  the  rate  of  decay,  at  least  for  the  two  cases 
shown,  was  the  same.  This  behavior  matches  that  of  a  driven-dissipative  system.  The  e- 
fold  relaxation  time  of  Eth,  te,  was  calculated  to  be  6.5  hours.  Burton  et  al.,  1975 
proposed  that  Dst  behaved  in  a  similar  way  and  could  be  described  by  a  simple 
differential  equation.  Burke  et  al.,  2009  applied  this  technique  to  modeling  Eth  using  £vs 
as  the  driver.  Since  Eth  is  related  to  the  exospheric  temperature  (Too)  linearly,  Too  can  also 
be  modeled  in  this  way.  The  following  sections  will  detail  the  Burke  et  al.,  2009  approach 
and  highlight  some  of  the  simplifying  assumptions  that  were  made  during  its 
development. 
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Figure  6:  Plots  of  magnetospheric  electric  field  8ys  (black)  and  the  natural 
logarithm  of  Eth  sw  (red)  for  the  disturbance  on  JD  204-211,  2004.  Vertical  lines 
mark  times  of  electric  field  decrease.  The  slanted  blue  lines  have  the  same  slopes 
indicating  that  Eth  sw  decays  exponentially  when  Sys  turns  off.  The  estimated  e-fold 
relaxation  time  is  6.5  hrs.  (Adapted  from  Burke  et  al.,  2009) 


Observed  Data 

Burke  et  al.,  2009  used  measured  thermospheric  orbit-averaged  density  and  orbit- 


averaged  altitude  from  the  GRACE  satellite  as  ground  truth  data.  These  data  were  used  to 
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calculate  exospheric  temperature  using  a  quadratic  fit  to  the  Jacchia  1977  model  (Burke, 
2008),  namely 


=Ysai(h)pi  M 

i= 0 


(6) 


where  Tm  is  the  exospheric  temperature  and  pl(/i)  is  the  orbit  averaged  neutral  density  in 
g/cm3  raised  to  the  ith  power.  The  term  aj(h)  is  a  coefficient  described  by  the  matrix 
equation 
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where  h  is  the  orbit  averaged  altitude  in  km.  Burke  (Burke,  et  al.,  2009)  took  an  orbit 
average  of  density  and  height  before  calculating  the  orbit  averaged  exospheric 


temperature. 

Once  the  exospheric  temperature  is  calculated,  the  total  energy  of  the 
thermosphere  can  be  calculated  using  the  empirical  formula 


Eth  ( h  >  100 km)  =  5.365  xlO17  +(8.727xl013)7;  (8) 
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Where  Eth  is  the  energy  of  the  thermosphere  in  Joules  and  Tm  is  the  orbit-averaged 
exospheric  temperature.  These  values  of  Eth  and  Tm  were  used  as  the  “observed”  data  for 
comparison  with  the  results  of  the  driven-dissipative  model  (Burke,  et  al.,  2009). 

Differential  Equations 

As  shown  in  Figure  6,  Eth  responds  to  changes  in  magnetospheric  electric  field  in 
a  way  reminiscent  of  a  driven-dissipative  thermodynamic  system.  Burton  et  al.,  1975, 
suggested  that  Dst  evolves  in  a  similar  way  and  developed  a  differential  equation  for  the 
pressure  corrected  Dst  (Dst*) 


dDst *  Dst *  (9) 

—  =  CCd£j 

dt  trc 

where  aD  is  the  coupling  coefficient,  s,  is  the  interplanetary  magnetic  field  magnitude 
and  trc  is  the  relaxation  time  constant  of  the  ring  current.  Dst*  is  defined  as  Dst*  = 

Dst  —  b^JPsw  +  c  where  b  and  c  are  empirical  constants  and  Psw  is  the  dynamic  pressure 
of  the  solar  wind  (Burton,  et  al.,  1975).  The  term  “driven-dissipative  system”  is 
illustrated  by  the  form  of  Equation  (9).  The  term  aDsI  models  the  driver  of  energy  input 

D* 

to  the  system  and  the  term  —  models  the  dissipation  of  energy  from  the  system. 

trc 

Burke  et  al.,  2009  used  this  approach  to  create  a  differential  equation  for  E*  using 
magnetospheric  electric  field  as  the  driver.  Using  this  equation  along  with  the  linear 
relationship  between  Eth  and  Too  from  Equation  (8)  yields  a  similar  equation  for  Too.  To 
simplify  the  model,  E*  and  Too  were  broken  into  two  independent  components,  one  due  to 
the  EUV  radiation  and  one  due  to  the  solar  wind  given  by 
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(10) 


E,i,  Eth  uv  +  E,i,  sw 


T  =T  +T 

Aco  ^  aaUV  x  SW 


(11) 


where  Eth  uv  =  6.1  x  10 17 J  and  Tcx>  uv  =  850  K,  both  considered  constant.  Then  Ethsw 
and  Too  sw  were  both  modeled  using  the  differential  equations 


dF 

SW 

dt 
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SW 
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where  aE  and  aT  are  the  coupling  constants  and  for  thermospheric  energy  and  exospheric 
temperature  respectively,  r  is  the  relaxation  time  constant,  the  same  for  both  parameters, 
and  sys  is  the  magnetospheric  electric  field  calculated  from  solar  wind  data. 

Equations  (12)  and  (13)  can  be  solved  numerically  for  any  time  in  the  future  using 
the  simple  Euler  method. 


E„„  (<„, )  =  E„sw  (t, )  +  SI  [  ccEsvs  (t, )  - 


(14) 


looSF  (^n+l)  ^oo  Sw{}n)Jr^- 


ai  F  s  (/« ) 


cSW  (/«  ) 


(15) 


Burke  et  al.,  2009  used  a  time  step  (At)  of  1  hour  in  their  analysis. 


32 


Relaxation  Constant 


The  relaxation  constant  (x)  is  defined  as  the  e-fold  relaxation  time  of  Eth  and 
determined  by  the  linear  fit  to  a  plot  of  the  natural  logarithm  of  Eth  sw  during  periods  of 
low  £[/5  after  storming  periods  seen  in  Figure  6.  Using  two  relaxation  periods  in  2004,  the 
constant’s  value  was  determined  to  be  r  «  6.5  hrs.  This  value  is  the  same  when  applied 
to  model  either  Eth  or  Too. 

Coupling  Constant 

Burke,  et  ah,  2009  used  comparisons  with  GRACE  data  from  JD  150-230,  2004 
to  determine  the  value  of  the  coupling  constant  for  thermospheric  energy,  aE. 


ar  =  5.5x10 


is 

hr*mV 


Using  this  value  good  agreement  was  shown  between  modeled  thermospheric  energy 
using  Equation  (14)  and  GRACE-derived  thermospheric  energy  using  Equation  (8) 
during  two  storming  periods  in  2004,  as  shown  in  Figure  7. 

Using  the  relationship  between  Eth  and  Tm  shown  in  Equation  (8),  the  coupling 
coefficient  for  exospheric  temperature  was  found  to  be: 


Using  this  value,  relatively  good  agreement  was  shown  between  modeled,  Equation  (15), 
and  GRACE-derived,  Equation  (6),  exospheric  temperatures  during  two  storming  periods 
in  2004,  as  shown  in  Figure  8. 
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Figure  7:  Comparison  of  £^5  (black),  modeled  Ethsw  (blue),  and  Ethsw  values 
inferred  from  GRACE  measurements  (red  dots)  plotted  as  functions  of  Universal 
Time  during  the  magnetically  disturbed  periods  of  July  (top)  and  November 
(bottom)  2004.  Ethsw  is  plotted  in  units  of  1016  J  (Burke,  et  al.,  2009). 


Figure  8:  Modeled  T ^  sw  (blue)  and  values  inferred  from  GRACE  measurements 
of  orbit  averaged  neutral  density  (red  dots),  plotted  as  functions  of  Universal  Time 
during  July  (top)  and  November  (bottom)  2004.  T^sw  was  approximated  by 
subtracting  850  K  from  GRACE-based  estimates  of  (Burke,  et  al.,  2009). 
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To  further  test  the  validity  of  the  coupling  constant  aE,  Burke  et  al.,  2009 
compared  the  term  aEevs,  which  represents  the  rate  at  which  energy  is  input  into  the 
thermosphere  from  the  solar  wind,  to  predictions  from  the  independent  W5  model 
(Wiemer,  2005).  The  W5  model  uses  IMF  and  solar  wind  measurements  to  predict  the 
Poynting  flux  into  the  ionosphere.  By  integrating  this  flux  over  the  polar  caps,  the  total 
rate  of  power  input  to  the  ionosphere  can  be  determined  and  compared  with  the 
predictions  from  the  term  aEsvs  in  the  driven-dissipative  model  (Burke,  et  al.,  2009). 
Figure  9  shows  that  during  two  storm  periods  in  2004  these  two  independent  models 
produce  similar  results,  validating  the  value  for  aE. 
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Figure  9:  Comparison  of  storm  time  power  into  the  global  thermosphere 
predicted  by  the  W5  model  (red)  and  aE £vs  (black)  plotted  as  functions  of  UT 
during  July  (top)  and  November  (bottom),  2004  (Burke  et  al.,  2009). 
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III.  Methodology 


Overview 

The  success  of  the  driven-dissipative  model  in  predicting  Toosw  and  Ethsw  for  the 
storm  periods  in  2004  is  promising  but  to  establish  general  applicability,  a  larger  sample 
of  storms  needs  to  be  studied.  This  section  outlines  the  methodology  used  to  test  the 
driven-dissipative  model  in  this  thesis.  First,  a  general  schematic  of  the  model  is 
discussed.  Then,  the  procedures  utilized  to  determine  the  observed  thermospheric 
temperatures  and  magnetospheric  electric  field  values  used  in  the  model  are  outlined. 
Second,  the  model’s  governing  equation  is  developed  and  solved.  Next,  the  procedure 
used  to  determine  optimal  coupling  and  relaxation  constant  values  for  each  storm  is 
presented.  Finally,  a  method  to  convert  model  temperature  values  to  model  density  values 
is  discussed. 

Model  Schematic 

A  general  schematic  of  the  model  is  shown  in  Figure  10.  The  model  uses 
observed  data  from  the  GRACE  satellite  to  derive  orbit-average  T\a  values  which  are 
used  as  the  data  the  model  attempts  to  replicate.  Observed  solar  wind  data  from  the  ACE 
satellite  is  used  to  calculate  magnetospheric  electric  field  magnitudes  which  serve  as  the 
driver  of  energy  input  to  the  thermosphere  in  the  model.  The  governing  equation  is  then 
solved  using  one  of  three  UV  methods  and  an  error  minimization  routine  which  selects 
values  for  the  coupling  and  relaxation  constants  for  each  storm.  This  process  results  in 
model  T1/2  data  for  each  method  and  each  storm.  Model  T1/2  data  is  then  converted  to 
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model  density  values  via  the  J77  model.  Each  step  is  explained  in  subsequent  sections  of 
this  chapter. 


Figure  10:  Schematic  of  the  Driven-Dissipative  model  used  in  this  Thesis.  Green  text  indicates 
observed  model  inputs.  Red  text  indicates  model  outputs. 
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Observed  Data 


Burke,  2011  conducted  a  study  of  38  geomagnetic  storms  from  2002  through 
2008  using  the  methods  outlined  above  from  Burke  et  al.,  2009.  The  storms  were  selected 
based  on  the  availability  of  the  solar  wind  data  necessary  to  compute  magnetospheric 
electric  field  values.  These  same  storm  periods  were  used  in  this  thesis,  with  one 
exception.  One  of  the  storms  used  by  Burke,  201 1  (Julian  Date  168-170,  2003)  did  not 
have  solar  wind  data  available  at  the  one-minute  time  cadence  used  for  this  thesis.  This 
storm  was  replaced  by  a  storm  from  Julian  Date  94  -  98,  2004  which  was  not  studied  by 
Burke,  2011.  The  storm  start  times,  end  times,  and  storm  types  (CME  or  CIR-driven)  for 
storms  used  in  this  thesis  are  listed  in  Table  1.  The  following  sections  outline  how  the 
observed  data  in  the  model  was  obtained. 

Storm  Period 

The  start  of  a  geomagnetic  storm  is  usually  defined  in  part  by  an  increase  in  solar 
wind  speed  and/or  density  coupled  with  a  southward  Z-component  of  the  solar  wind 
magnetic  field  (Bz  south).  Using  the  initial  days  listed  in  Burke,  2011  for  each  storm 
period  as  a  starting  point  the  time  of  this  increase  was  determined  for  each  storm.  The 
storm  starting  time  was  defined  as  the  last  time  the  GRACE  satellite  crossed  the  equator 
on  an  ascending  pass  prior  to  the  increase  in  solar  wind  speed  and/or  density.  The  end  of 
the  storm  period  was  generally  defined  as  the  final  day  listed  by  Burke,  201 1.  In  some 
cases,  that  time  was  clearly  after  both  the  magnetospheric  electric  field  and  GRACE- 
derived  exospheric  temperature  had  recovered  to  a  state  of  quasi-equilibrium  near  pre¬ 
storm  levels.  In  these  cases  the  storm  end  time  was  adjusted  backwards  to  match 
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Table  1:  Storm  Periods  Studied 


Year 

Day 

Storm  Start 

Hour  Min 

Sec 

Day 

Storm  End 

Hour  Min 

Sec 

Storm 

Type 

246 

18 

16 

59 

248 

0 

0 

0 

CME 

2002 

250 

13 

35 

35 

252 

0 

0 

0 

CME 

2002 

272 

23 

56 

57 

275 

16 

48 

0 

CIR 

2002 

276 

9 

47 

28 

278 

0 

0 

0 

CIR 

2002 

296 

22 

24 

0 

299 

0 

0 

0 

CIR 

2002 

324 

14 

23 

18 

327 

0 

0 

0 

CIR 

2003 

149 

12 

22 

43 

151 

0 

0 

0 

CME 

2003 

229 

14 

1 

55 

232 

0 

0 

0 

CIR 

2003 

324 

6 

54 

25 

326 

0 

0 

0 

CME 

2004 

22 

0 

37 

57 

24 

0 

0 

0 

CME 

2004 

94 

0 

54 

55 

98 

0 

0 

0 

CIR 

2004 

204 

9 

53 

59 

210 

0 

0 

0 

CME 

2004 

208 

22 

10 

40 

210 

0 

0 

0 

CME 

2004 

243 

1 

12 

42 

246 

0 

0 

0 

CIR 

2004 

312 

9 

44 

47 

314 

0 

0 

0 

CME 

2004 

314 

10 

22 

30 

316 

0 

0 

0 

CME 

2005 

127 

17 

35 

50 

130 

0 

0 

0 

CIR 

2005 

135 

4 

20 

33 

136 

0 

0 

0 

CME 

2005 

148 

4 

20 

25 

152 

0 

0 

0 

CME 

2005 

163 

8 

4 

12 

165 

0 

0 

0 

CIR 

2005 

236 

0 

53 

0 

239 

0 

0 

0 

CME 

2005 

243 

8 

26 

32 

245 

0 

0 

0 

CIR 

2006 

98 

10 

5 

11 

100 

19 

33 

36 

CIR 

2006 

103 

4 

29 

30 

107 

0 

0 

0 

CIR 

2006 

348 

14 

1 

20 

350 

0 

0 

0 

CME 

2007 

142 

7 

30 

53 

148 

9 

36 

0 

CIR 

2007 

191 

20 

17 

30 

193 

0 

0 

0 

CIR 

2007 

195 

7 

18 

0 

197 

0 

0 

0 

CIR 

2007 

218 

12 

52 

5 

221 

0 

0 

0 

CIR 

2007 

298 

9 

50 

7 

305 

0 

0 

0 

CIR 

2007 

323 

17 

29 

20 

325 

21 

14 

24 

CIR 

2007 

351 

3 

12 

50 

356 

0 

0 

0 

CIR 

2008 

31 

13 

9 

13 

37 

0 

0 

0 

CIR 

2008 

68 

9 

55 

30 

71 

0 

0 

0 

CIR 

2008 

86 

2 

17 

0 

91 

0 

0 

0 

CIR 

2008 

166 

17 

38 

35 

170 

0 

0 

0 

CIR 

2008 

194 

0 

8 

55 

195 

0 

0 

0 

CIR 

2008 

247 

2 

40 

2 

250 

0 

0 

0 

CIR 
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observations.  Temperatures  derived  from  observed  GRACE  density  and  height  values  are 
referred  to  as  observed  temperature  values  for  the  remainder  of  the  thesis.  Figure  1 1 
shows  an  example  of  a  typical  storm  period. 
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Figure  11:  Example  of  a  storm  period,  defined  as  the  start  time  of  the  GRACE 
orbit  just  prior  to  the  initial  electric  field  rise  until  both  the  electric  field  and  the 
observed  Ti/2  have  returned  to  quasi-equilibrium  near  their  pre-storm  values. 
Observed  GRACE  T1/2  (top)  and  magnetospheric  electric  field  data  (bottom)  are 
shown  as  functions  of  Julian  date  (JD),  2007  where  JD  is  counted  from  1  Jan, 
2007.  The  red  vertical  lines  indicate  the  storm  start  time  (left)  and  storm  end  time 
(right). 


Exospheric  Temperature 

Once  the  storm  periods  were  determined,  data  from  the  GRACE  A  satellite  was 
used  to  calculate  the  exospheric  temperature,  Too.  The  GRACE  A  satellite  was  in  a  polar 
orbit  at  altitudes  from  455-534  km  during  the  period  from  2002-2008.  The  GRACE  data 
set  used  for  this  thesis  was  calibrated  by  Sutton,  201 1  with  thermospheric  parameters 
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averaged  into  3  degree  latitude  bins.  For  this  research  neutral  density,  altitude,  latitude 
and  local  time  were  used.  Burke  et  al.,  2009  used  an  earlier  calibrated  version  of 
GRACE  data  (Burke,  2011)  that  was  not  averaged  into  latitude  bins.  His  approach  was  to 
orbit-average  the  density  and  altitude  from  raw  GRACE  data  first,  and  then  calculate  the 
orbit-averaged  T*  using  Equations  (6)  and  (7).  Wise  et  al.,  2012  showed  that  this  method 
of  orbit  averaging  produced  inaccurate  results.  The  GRACE  satellite  orbit  is  slightly  non¬ 
circular  with  an  apogee  about  20km  higher  than  perigee.  Orbital  dynamics  dictate  that  the 
satellite  moves  slower  near  apogee  than  perigee.  This  results  in  underestimates  of  orbit- 
average  density  and  orbit-average  heights  higher  than  the  time-independent  average. 
Combined,  these  factors  lead  to  underestimates  of  orbit-averaged  Tx  when  it  is  calculated 
from  orbit-averaged  heights  and  densities  (Wise,  et  al.,  2012). 

Wise  et  al.,  2012  showed  that  a  more  physically  accurate  way  to  calculate  orbit- 
averaged  Too  is  to  calculate  Tx  in  each  GRACE  3 -degree  latitude  bin  and  then  average  the 
result.  In  this  thesis  I  used  the  approach  of  Wise,  2012  and  computed  exospheric 
temperature  for  each  latitude  bin  prior  to  the  orbit  averaging,  resulting  in  an  observed  Tx 
for  each  GRACE  latitude  bin.  The  change  in  technique  does  impact  the  resulting 
observed  orbit-average  temperature  values  and  results  are  shown  in  section  IV. 

The  method  to  compute  exospheric  temperature  used  by  Burke  et  al.,  2009 
(Equations  (6)  and  (7) )  is  taken  from  Burke’s  (2008)  quadratic  fit  to  the  J77  model.  This 
fit  was  developed  to  be  accurate  only  within  the  ranges  of  300  -  500  km  in  altitude  and 
700  -  2000  K  in  Tx  which  leaves  open  the  possibility  that  the  fit  is  not  sufficiently 
accurate  over  the  entire  temperature  range  present  in  the  38-storm  sample  listed  in  Table 
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1.  To  quantify  and  correct  this  possible  source  of  error  an  interpolation/iteration 
technique  was  developed  to  produce  an  exospheric  temperature  that,  coupled  with  the 
GRACE  altitude,  produces  the  observed  GRACE  density  via  the  J77  model.  To  begin, 
tables  of  data  from  the  J77  model  were  generated  using  a  Fortran  code  written  by  David 
Huestis  in  1999  and  provided  by  John  Wise  at  the  Air  Force  Research  Laboratory.  Data 
tables  list  densities  as  a  function  of  altitude  for  a  given  exospheric  temperature.  Tables 
were  generated  for  exospheric  temperatures  from  500-2000K  with  a  resolution  of  100K, 
listing  densities  for  altitudes  of  300km- 1000km  with  a  resolution  of  1km. 

Data  from  the  J77  tables  were  used  to  create  a  3-D  grid  of  data  giving  density  for 
a  specified  Too  -  altitude  pair.  The  temperature  and  altitude  ranges  chosen  ensure  that  all 
of  the  observed  GRACE  data  fit  inside  the  data  grid.  With  the  data  grid  as  a  basis,  a 
density  can  be  generated  using  any  specified  Too,  altitude  pair  by  interpolating  between 
the  data  points.  For  this  thesis,  cubic  spline  interpolation  (Press,  et  al.,  2007)  was  used  via 
MATLAB’s  interp2  function.  To  generate  observed  exospheric  temperatures  from 
observed  GRACE  heights  and  densities  an  iterative  technique,  the  Nelder-Mead  simplex 
direct-search  method  (Lagarias,  et  al.,  1998),  was  used.  Starting  at  an  initial  guess  for  Too, 
here  800K,  the  search  method  iterates  over  T®  values  until  a  T*,  is  found  that  minimizes 
the  relative  error  (to  a  tolerance  of  10~4  %)  between  the  observed  density  and  the 
interpolated  density  when  paired  with  the  observed  altitude.  The  Nelder-Mead  method  is 
detailed  in  Appendix  B. 

Thermospheric  models  like  J70,  J77  and  HASDM  use  a  global  temperature 
parameter  to  account  for  the  EUV  contribution  to  the  thermosphere’s  energy  budget.  The 
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J77  uses  the  arithmetic  mean  temperature,  T\a,  defined  as  the  average  of  the  daytime 
maximum  T«,  and  the  nighttime  minimum  Tx  at  a  given  time.  While  Too  characterizes  the 
thermosphere  at  a  specific  location  and  time,  T  m  is  general  parameter  that  removes  the 
diurnal  variation  of  Too  and  characterizes  the  state  of  the  thermosphere  as  a  whole.  The 
formula  from  the  J77  model  (Jacchia,  1977)  shown  in  Equation  (18)  was  used  to  convert 
the  observed  Too  values  from  the  GRACE  data  into  T\a. 


T  T,  (18) 

1/2  D(S,  (/),H) 

The  conversion  factor  D  is  a  function  of  solar  declination  angle  ( 8 ),  latitude  (</>)  and  solar 
hour  angle  (77)  and  given  by  Equation  (19). 


<7  . 


D  =  l  +  Cl-Sin(^)  +  c2|/(/7)  ^ 


1 


cos 


w 


(19) 


where: 


m= 


cos 


;{«+/>) 


+  c, 


cos[3(77  +  ^)  +  ^] 


Cj  =0.15,  c2  =  0.24,  c3  =0.08,  y9  =  -60°,  e  =  23.44°,  j  =  -75° 


Latitude  was  taken  as  the  mean  location  of  the  GRACE  satellite  in  each  latitude  bin.  The 
hour  angle  H  is  simply  the  mean  local  time  of  each  latitude  bin,  converted  to  an  angle 
counted  from  local  noon  via  the  formula: 
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H  =  [Local  Time(hours)  -12)xl5 


(20) 


The  solar  declination  angle  (S)  was  calculated  using  the  time  of  each  data  point  via  the 
method  described  by  Meeus,  1991  outlined  in  Appendix  A. 

Once  observed  Too  values  were  converted  to  T\a  via  Equation  (18)  in  each  latitude 
bin,  the  observed  Tm  data  was  orbit  averaged.  The  start  of  each  orbit  was  defined  as  the 
equator  on  each  ascending  pass  of  the  GRACE  satellite.  The  end  of  each  orbit  was 
defined  as  the  point  just  before  the  equator  on  each  ascending  pass.  All  of  the  Tm  values 
for  each  orbit  were  averaged  to  produce  a  single  value  for  each  orbit,  and  the  time  for 
each  orbit-averaged  value  was  defined  as  the  time  of  the  start  of  the  orbit.  The  resulting 
orbit-averaged  T|/2  values  and  times  were  used  as  the  observed  data  the  model  attempts  to 
replicate. 

Magnetospheric  Electric  Field 

The  main  source  of  energy  for  the  thermosphere  during  geomagnetic  storms  is  the 
solar  wind  which  couples  to  the  thermosphere  via  the  magnetospheric  electric  field,  sys- 
Using  solar  wind  data  from  the  Advanced  Composition  Explorer  (ACE)  satellite,  sys  can 
be  calculated  in  near  real-time  using  a  version  of  the  Volland-Stem  model  originally 
formulated  by  Ejiri,  1978  and  modified  by  Burke,  2007.  The  ACE  satellite  is  located  at 
the  LI  Lagrange  point  between  the  Sun  and  the  Earth  which  is  roughly  one  hour 
upstream  of  the  Earth  in  the  solar  wind  flow.  For  this  thesis  ACE  data  that  had  already 
been  time  shifted,  meaning  the  time  stamp  on  the  data  was  adjusted  by  roughly  one  hour 
to  account  for  the  transit  time  to  earth,  was  utilized.  The  exact  amount  of  time  adjustment 
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depends  on  the  current  solar  wind  velocity.  This  data  is  available  at  a  one-minute  time 
cadence  from  the  NASA  OMNIWeb.  Solar  wind  pressure,  Pswi  and  velocity,  Vsw,  were 
obtained  from  the  ACE  Solar  Wind  Electron  Proton  Alpha  Monitor  (SWEPAM)  while 
the  Y  (By)  and  Z  (Bz)  components  of  the  solar  wind  magnetic  field  were  obtained  from 
the  ACE  Magnetic  Field  Experiment  (MFE)  sensor.  All  calculations  and  data  use  GSM 
coordinates,  illustrated  in  Figure  2. 

Using  Burke’s  (Burke,  2007)  formulation,  the  magnetospheric  electric  field 
magnitude  can  be  calculated  using  the  relation 

F  _  ^  pc  (21) 

2  LyRe 

The  denominator  in  Equation  (21)  gives  the  width  of  the  magnetosphere  in  the  Y 
direction.  RE  is  the  radius  of  earth  and  LY  is  the  distance  to  the  magnetopause  in  the  Y 
direction,  in  earth  radii,  calculated  using  the  solar  wind  pressure  with  Equation  (22). 


14.4 

yjpsw  (nPa) 


(22) 


The  numerator  in  Equation  (21),  @PC,  is  the  cross-polar  cap  electric  potential. 
Siscoe  et  al.,  2002  built  on  the  Hill  model  (Hill,  1984)  and  developed  a  formula  for  Opc 
using  the  magnetospheric  saturation  potential,  @s,  and  the  magnetospheric  convection 
potential,  Oe  as  inputs. 
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(23) 


oE+os 

Equation  (23)  shows  that  &s  serves  as  a  limiting  value  for  0PC.  In  other  words,  when 
0E  «  0PE  ^  ^  while  when  ^  ^  0pc  ^  ^5  (Hill,  et  al.,  1976). 

@s  is  the  potential  that  drives  region  one  currents  in  the  magnetosphere  which 
create  magnetic  fields  that  weaken  the  earth’s  magnetic  field  at  the  magnetopause 
(Siscoe,  et  ah,  2002).  It  can  be  calculated  using  the  solar  wind  dynamic  pressure  (Psw) 
using  the  formula 


1600 (nPa)  (24) 

s  Y 

^ p 

where  Ip  is  the  effective  Pedersen  conductance  of  the  polar  cap,  here  approximated  as  a 
constant  27p  =  10  mho  (Burke,  2007) .  We  see  that  increased  solar  wind  pressure  results 
in  a  greater  magnetospheric  saturation  potential. 

0E  results  from  magnetic  reconnection  processes  at  the  magnetopause 
(Boudouridis,  et  ah,  2004).  It  can  be  calculated  using  the  solar  wind  velocity  and 
magnetic  field  data  via  Equation  (25) 


0O  +  LgVswBt  sin 


(25) 


where  the  first  term,  &0,  is  a  residual  potential  due  to  viscosity  in  the  low-latitude 
boundary  layer  (Burke,  2007),  (Kennel,  1995).  Burke  found  that  &0  typically  ranges 
between  20  and  30  kV  (Burke,  et  ah,  1999)  and  for  this  research  the  value  of  @0  was  set 
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Q 

at  25  kV.  The  second  term  in  Equation  (25),  VSWBT  sin2  -,  gives  the  magnitude  of  the 
interplanetary  electric  field  as  developed  by  Sonnerup,  1976,  where  Vsw  is  the  solar 
wind  velocity,  BT  =  J By  +  B2  and  0  =  cos-1  — ,  the  interplanetary  electric  field  clock 

Bj1 

angle  in  the  Y-Z  plane.  The  interplanetary  electric  field  multiplied  by  Lg,  the  width  of  the 
space  (in  Earth  radii)  through  which  geoeffective  solar  wind  streamlines  must  pass  to 
reach  the  dayside  magnetopause  (Burke,  2007),  gives  the  interplanetary  electric  potential. 
Typically  LG  values  between  3-4  Earth  radii  and  in  this  research  the  approximation  is 
made  that  LG  =  3.5,  a  constant,  as  suggested  by  Burke  et  al.,  1999. 

Data  from  the  ACE  satellite  is  occasionally  either  bad  or  missing.  When  missing 
or  bad  data  made  reliable  electric  field  values  impossible  to  calculate  directly, 
interpolation  was  used  between  the  nearest  good  data  points  to  fill  in  the  gap  and  ensure 
good  electric  field  values  existed  for  each  minute  during  storm  time.  Because  storm 
periods  were  selected  based  on  relatively  good  ACE  data  availability,  the  amount  of 
interpolation  was  kept  to  a  minimum.  None  of  the  storms  studied  had  contiguous  gaps  in 
ACE  data  of  longer  than  four  hours  and  none  of  the  storms  studied  had  missing  ACE  data 
at  the  time  of  peak  observed  temperatures.  Figure  12  shows  an  example  of  the  ACE  solar 
wind  data  and  the  magnetospheric  electric  field  magnitude  calculated  using  the  above 
formulation  for  a  geomagnetic  storm  in  July,  2004. 
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Figure  12:  Illustration  of  ACE  solar  wind  data  and  the  resultant  magnetospheric 
electric  field  magnitude.  From  top  to  bottom  the  figure  shows  the  y  and  z 
components  of  the  interplanetary  magnetic  field  (By  and  Bz)  in  nano-Teslas,  the 
solar  wind  pressure  (P)  in  nano-Pascals,  the  solar  wind  velocity  (V)  in  km/s,  and 
the  magnetospheric  electric  field  magnitude  (E  field)  in  miliVolts/meter  as 
functions  of  modified  Julian  date  (JD),  2004  where  JD  is  counted  from  1  Jan, 
2004. 


Governing  Equation 

Burke’s  original  model  (Burke,  et  al.,  2009)  assumed  that  the  UV  contribution  to 
the  exospheric  temperature  was  constant  throughout  the  storm  period.  It  is  true  that 
during  storm  time,  the  geomagnetic  contribution  to  thermospheric  energy  and  therefore 
exospheric  temperature  is  much  more  variable  than  the  UV  contribution.  However,  since 
the  sun’s  UV  irradiance  can  change  on  short  time  scales  as  well,  a  more  realistic  model 
would  allow  the  UV  contribution  to  vary  during  storm  time.  Since  the  J77  model 
accounts  for  the  UV  contribution  via  the  arithmetic  mean  temperature  (T1/2),  T1/2  is 
modeled  via  a  new  differential  equation  governing  its  time  rate  of  change. 
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Since  T\a  differs  from  Too  by  only  a  conversion  factor  shown  in  Equation  (18), 
T 1/2  can  be  expressed  as  the  sum  of  UV  and  solar  wind  contributions  just  as  Too  was  by 
Burke  in  Equation  (11). 


T 

1/2 


—  T  TTV  +  T 


uv 


,  sir 


(26) 


Likewise  the  time  rate  of  change  of  T\n  is  simply  the  sum  of  the  time  rates  of  change  of 
its  components. 


dT 
_ _ 1 /_ 

dt 


^1/2 uv_  ,  dTyisw 

dt  dt 


(27) 


Just  prior  to  storm  time  the  thermosphere  is  taken  to  be  at  equilibrium  with 


dTi/2 

dt 


=  0  and 


T 

1/2 


_  rji  0  _  rj-i  0 


UV 


+r 


sw 


(28) 


Equation  (28)  expresses  the  pre-storm  equilibrium  arithmetic  mean  exospheric 
temperature  ( T^2 )  as  the  sum  of  the  equilibrium  UV  and  solar  wind  contributions.  In 
general,  outside  geomagnetic  storming  periods,  the  UV  contribution  to  thermospheric 
energy  is  much  larger  than  the  solar  wind  contribution  (Figure  1),  which  indicates  that 
ti/2uv  »  ti/2sw  suggesting  the  approximation  T?/wv  =  T?/2. 

Using  Burke’s  (Burke,  et  al.,  2009)  expression  for  given  by  Equation  (13) 

and  substituting  d'lj2sw  =  Tj/2  —  7i°/2,  Equation  (27)  becomes 


dT 

_ 1A 

dt 


dT  uv  T  (t)-T° 

+  asvs  ( t )  — ^2 — 

dt  VS  '  ’  T 


(29) 
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where  a  is  a  coupling  constant  linking  the  magnetospheric  electric  field  (svs)  to  T \/2  and 
x  is  a  relaxation  constant.  Solving  Equation  (29)  using  Euler’s  method  provides  a  model 
of  T 1/2  as  a  function  of  time  to  compare  with  the  observed  data  from  GRACE. 


UV  contribution  to  T1/2 

Equation  (29)  shows  that  any  solution  for  T1/2  depends  on  the  ways  in  which  both 
Ti/iuv  and  7i°/2  are  treated.  They  can  be  treated  as  constants  through  the  storm  period  or 
allowed  to  vary.  If  they  are  variable,  a  method  of  calculating  their  value  must  be  selected. 
The  following  sections  outline  the  three  ways  in  which  T^uy  and  Ty2  are  treated  for 
this  research. 

Method  One 

The  simplest  way  to  treat  T1/2uv  and  T®/2  is  to  approximate  them  as  constants 
through  the  storm  period.  This  is  the  method  Burke  et  al.,  2009  used  in  their  original 
model.  In  method  one,  dTl/2t/t'  =  0  at  all  times  and  Equation  (29)  becomes 


dT  .  .  T  (t)-T 

-7r(t)  =  aers(‘)-  —  ' 

Clt  T 


(30) 


Solving  Equation  (30)  via  the  Euler  method  results  in  a  time  dependent  formula  for  ^1/2 
shown  in  Equation  (31). 


T 

1/2 


(t  +  At)-T/2(t)  +  At 


as, 


vs 


(0 


T 

_ 1/2 


(')-/: 


T 


(31) 
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In  method  one,  TyZ  was  defined  as  the  mean  of  the  observed  arithmetic  mean  exospheric 
temperature  from  the  8  GRACE  orbits  (=12  hrs)  prior  to  the  storm  start  time  and 
considered  constant  throughout  the  storm.  A  time  step  (At)  of  1  minute  was  used  to  match 
the  cadence  of  the  ACE-derived  magnetospheric  electric  field  data. 

Method  Two 

In  method  two  the  pre-storm  equilibrium  temperature,  T®/2,  was  still  considered  to 
be  constant  throughout  the  storm  period  and  defined  identically  to  method  one.  However, 
was  a"owed  to  be  non-zero  in  Equation  (29).  Solving  Equation  (29)  with  a  non¬ 
zero  dTl/2W  via  the  Euler  method  with  a  one-minute  time  step  results  in  the  time- 

dt  ^ 

dependent  expression  for  Tl/2  used  in  method  two. 


T  (f  +  A/^ —  T  (y^  +  At 


dT 


uv 


dt 


(0  +  ccTsvs  (V) 


T  (t)-T° 

1/2  V  /  1/2 


(32) 


To  account  for  the  variation  in  the  UV  contribution  to  the  exospheric  temperature, 
the  J77  model  calculates  T1/2  uv  as  a  function  of  the  F10.7  index  using  the  formula 


TrW  (0  =  5.48  (/•„„ )“  + 101  ,&(Fm f  (33) 

where  F10.7  is  simply  the  daily  value  of  the  F10.7  index  and  F  10.7a  is  an  162-day  averaged 
value  of  the  F10.7  index.  Using  the  results  of  Equation  (33),  the  time  rate  of  change  of 
T1/2uv  was  calculated  for  each  minute  during  the  storm  period  using  Equation  (34). 
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(34) 


dTin uv  /  x  T  uv(next  observation)  -  T  uv  (prior  observation) 
dt  24  hrs 

Equation  (34)  results  in  a  value  for  for  each  minute  during  the  storm  period  with 

units  of  The  J77  formulation  for  T^uv  in  Equation  (33)  results  in  values  for 

that  are  constant,  but  not  necessarily  zero,  for  24  hour  periods  between  F10.7  observations 
at  20Z  each  day. 

Method  Three 

In  method  three,  was  allowed  to  vary  using  Equations  (33)  and  (34)  in  the 

same  way  as  method  two.  T®/2  was  also  allowed  to  vary  by  approximating 
T°i 2  =  T-[°2UV,  where  Ty2UV  is  a  corrected  version  of  T1/wv.  To  obtain  values  for 
Tyl uv(t)  at  a  one-minute  time  cadence,  the  ill  formula  for  Txj2UV  was  used  (Equation 
(33))  to  calculate  a  value  for  T1/2UV  at  20Z  each  day  and  then  interpolated  to  produce  a 
value  at  each  minute  during  the  storm  period.  Because  Equation  (33)  is  a  modeled  input, 
it  does  not  always  match  the  observed  value  of  T°/2,  as  defined  in  method  one,  at  the 
beginning  of  the  storm  period.  To  remove  this  discrepancy  a  correction  factor,  K,  was 
added  to  the  modeled  value  oiT1/2UV  at  each  time 

TZ,  (')=  r„ K  (35) 

where  T'i/2uv(t)  is  the  result  of  Equation  (33)  after  interpolation  and  K  is  given  by 
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(36) 


0 ,  observed  0 ,  modeled 

K  =  TUV  -T 


1/2 


1/2 


uv 


In  Equation  (36)  T^’°2l^yrveA\&  the  observed  pre-storm  equilibrium  temperature  defined  in 
method  one  and  T^’^ySled  is  the  modeled  value  of  the  exospheric  temperature,  using 
Equation  (33)  and  interpolation,  at  the  start  of  the  storm  period. 

The  approximation  T°/2  =  Tf°2UV  results  in  a  modified  version  of  Equation  (29) 
to  be  used  for  method  three,  given  by  Equation  (37). 
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(37) 


Solving  Equation  (37)  with  Euler’s  method  using  a  one-minute  time  step  results  in  the 
time-dependent  expression  for  ^1/2  used  in  method  three. 
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1/2 
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uv 
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(38) 


Orbit  Averages 

Because  the  observed  T  i/2  data  from  GRACE  that  the  model  is  attempting  to 
replicate  is  averaged  over  the  period  of  an  orbit,  the  modeled  data  needs  to  be  averaged 
over  the  same  time  periods  to  facilitate  direct  comparison.  After  modeled  T  i  /2  values  are 
calculated  for  each  minute  of  the  storm  period  using  one  of  the  three  methods  above  the 
modeled  values  were  averaged  over  the  same  time  periods  defined  earlier  by  the  GRACE 
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orbits.  The  time  stamp  for  each  averaged  value  is  taken  to  be  the  start  of  the  average 
period,  which  is  equivalent  to  the  start  of  the  GRACE  orbit  used  for  the  observed  data. 

Coupling  Constant  and  Relaxation  Constant 

The  last  piece  of  the  model  that  must  be  determined  is  the  value  of  the  coupling 

constant,  a,  and  the  relaxation  constant,  x,  in  Equations  (31),  (32),  and  (38)  for  methods 
one,  two  and  three  respectively.  The  values  of  both  a  and  x  were  considered  constant 
through  each  storm  period  but  were  allowed  to  have  different  values  for  each  storm 
period  and  for  each  method.  To  determine  the  optimal  value  of  a  and  x  for  each  storm 
period,  the  MATLAB  fminsearch  function  was  used  to  minimize  the  relative  root-mean- 
squared  (RMS)  error,  defined  by  Equation  (39),  between  the  observed  and  modeled 
values  of  T1/2  by  adjusting  the  values  of  a  and  x  . 


Relative  T,n  RMS  Error  = 

) 

Here  j°^erved  denotes  the  orbit-averaged  observed  T y2  values  derived  from  GRACE 
data,  T^°del  denotes  the  orbit-averaged  model  T1/2  values  using  one  of  the  methods 
described  above  and  N  is  the  number  of  data  points  during  the  storm  period.  The 
MATLAB  fminsearch  function  uses  the  Nelder-Mead  simplex  direct  search  method 
(Lagarias,  et  al.,  1998)  to  minimize  a  given  function.  The  algorithm  is  outlined  in 
Appendix  B.  Using  this  procedure  optimal  values  for  a  and  x  were  determined  for  each 
storm  period  and  each  method,  along  with  the  resulting  relative  RMS  error. 
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Density  Conversion 

The  main  purpose  of  the  driven-dissipative  model  is  to  show  the  relevance  of 
using  the  magnetospheric  electric  field  to  model  the  energy  input  to  the  thermosphere 
during  geomagnetic  storming.  This  is  accomplished  by  modeling  orbit-averaged  T  i/2. 
However,  since  there  are  no  published  temperature  errors  for  current  thermospheric 
models  such  as  HASDM,  the  temperature  errors  from  the  driven-dissipative  model  cannot 
be  directly  compared  to  existing  thermospheric  models.  In  order  to  facilitate  comparisons 
with  published  mean  HASDM  density  errors  of  17%  during  geomagnetic  storming 
(Marcos  et  al.,  2010)  the  model  T|/2  results  must  be  used  to  generate  model  density 
values.  In  an  operational  application,  this  would  be  done  by  using  the  driven-dissipative 
method  or  model  T|/2  output  in  a  current  thermospheric  model  such  as  HASDM  or  JB08. 
For  the  purpose  of  assessing  the  relationship  between  Ti/2  errors  and  density  errors,  the 
J77  model  can  be  used  to  generate  model  densities  from  model  Ty2  values. 

The  model  used  in  this  thesis  is  designed  to  minimize  the  error  between  observed 
and  modeled  orbit-average  Tj/2  values,  not  the  error  between  modeled  and  observed  T  i/2 
at  any  given  point  and  time.  Observed  T1/2  values  reflect  the  GRACE  satellite’s  latitude 
while  model  T1/2  values  do  not  exhibit  this  same  variation.  Instead  they  are  generated  by 
the  Equation  (38)  which  produces  T i/2  as  a  function  of  time  that,  while  not  matching  the 
latitudinal  variation  of  the  observed  T  i/2>  results  in  orbit-average  values  very  close  to 
those  observed.  This  relationship  is  shown  in  Figure  13.  The  observed  Ti/2  values  (blue 
line)  vary  as  a  function  of  latitude  (bottom  plot)  while  the  modeled  Tj/2  values  (red  line) 
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do  not.  Despite  this  difference,  orbit-averaged  observed  T\a  values  (blue  x’s)  agree  well 
with  model  orbit-averaged  T 1/2  values  (red  dots). 


Figure  13:  Illustration  of  the  latitudinal  dependence  of  observed  (blue)  and 
modeled  (red)  T1/2  values.  The  top  plot  shows  observed  T1/2  values  in  each  GRACE 
latitude  bin  (blue  line)  and  raw  modeled  T1/2  values  at  a  one  minute  time  cadence 
(red  line).  Orbit-averaged  observed  T1/2  values  are  shown  as  blue  x’s  and  orbit- 
averaged  model  T1/2  values  are  shown  as  red  dots.  The  bottom  plot  shows  that 
Latitude  of  the  GRACE  satellite  as  a  function  of  Julian  Date,  2004. 


Since  the  model  is  not  formulated  to  accurately  model  Ti/2  in  each  GRACE 
latitude  bin,  the  model  cannot  be  expected  to  accurately  model  density  in  each  GRACE 
latitude  bin.  This  means  that  the  orbit-averaging  method  used  previously,  where  values 
were  computed  for  each  latitude  bin  and  then  averaged  over  an  orbit,  cannot  be  used 
when  converting  model  Ti/2  into  model  density  values.  Instead,  to  generate  orbit-average 
model  densities,  orbit-average  model  T i/2  values  are  used.  Since  the  J77  model  tables  list 
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density  as  a  function  of  Too  and  height,  model  orbit-averaged  T1/2  is  converted  to  model 


orbit-averaged  Too  via  Equation  (40) 


^a>, mod  —  ^1/2, mod  '  ^  (40) 

where  fmmod  is  the  modeled  orbit-averaged  Too,  fy2 ,mod  is  the  modeled  orbit-averaged 
T 1/2,  and  D  is  the  observed  orbit-averaged  conversion  factor  from  the  J77  model 
calculated  via  Equation  (19). 

Once  Too, mod  was  calculated  it  was  paired  with  the  corresponding  observed 
GRACE  orbit-averaged  height  to  generate  model  orbit-averaged  density  via  interpolation 
within  the  ill  model  tables.  Using  this  process,  model  orbit-averaged  density  was 
calculated  for  each  storm  and  then  compared  to  the  observed  GRACE  orbit-averaged 
density  to  calculate  the  relative  density  RMS  error  via  Equation  (41) 


Relative  Density  RMS  Error  = 


(  P oh  P mod  ) 

Pob 


z 


(41) 


where  pob  is  the  observed  orbit-averaged  density,  pmod  is  the  modeled  orbit-averaged 
density  and  N  is  the  number  of  data  points. 
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IV.  Analysis  and  Results 


Chapter  Overview 

Section  IV  begins  with  a  comparison  of  the  method  used  by  Burke  et  al.,  2009  to 
derive  observed  exospheric  temperatures  from  GRACE  data  with  the  method  developed 
for  this  thesis.  Second,  results  are  presented  using  the  three  model  methods  for  the  38 
storms  in  the  sample.  Results  from  three  individual  storms  are  presented  in  detail  and  the 
results  of  the  different  methods  are  compared.  Model  results  are  compared  with  the 
results  from  Burke,  201 1  and  differences  are  discussed.  Next,  values  for  a  and  x 
determined  for  each  storm  are  fit  as  functions  of  Fio.7a  in  an  effort  to  make  the  model 
operationally  useful.  Results  of  the  fits  are  presented.  To  test  general  applicability,  the 
model  is  applied  to  two  storms  outside  the  original  sample  of  38.  Finally,  model 
temperature  values  are  converted  to  model  density  values  and  the  resulting  errors  are 
compared  with  published  HASDM  density  errors. 

Observed  Data 

Before  running  the  model,  observed  temperature  values  need  to  be  determined 
from  GRACE  data.  As  discussed  in  section  III,  Burke’s  approach  (2009)  was  to  calculate 
orbit-averaged  values  of  height  and  density  from  raw  GRACE  data  and  then  calculate  an 
orbit-averaged  T:x  value.  Wise  et  al.,  2012  showed  that  calculating  orbit-averaged  Tx 
from  orbit-averaged  density  and  height  data  produces  an  inaccurate  result  due  to  the  fact 
that  GRACE  orbits  are  slightly  non-circular.  Instead,  the  more  accurate  approach  is  to 
average  GRACE  density  and  altitude  into  3-degree  latitude  bins,  calculate  TX;  for  each 
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GRACE  latitude  bin,  and  then  orbit-average  the  resulting  values.  Wise’s  approach  to 
orbit-averaging  values  is  applied  here. 

After  calculating  orbit-averaged  density  and  height  values,  Burke  et  al.,  2009, 
applied  the  quadratic  fit  to  the  J77  model  given  in  equations  (6)  and  (7)  to  calculate  the 
orbit-averaged  T,x  values.  This  fit  was  developed  to  provide  good  results  for  heights 
within  the  range  from  300  to  500  km  and  for  Tx  from  700  to  2000  K  (Burke,  2008).  For 
the  storms  sampled  in  this  thesis,  GRACE  heights  range  from  455-534  km  and  Tx  ranged 
from  roughly  500K  -  1400K.  Since  the  GRACE  data  for  the  storm  sample  does  not  fit 
entirely  within  the  range  treated  well  by  the  quadratic  fit  a  test  was  run  to  determine  if  the 
quadratic  fit  would  produce  accurate  results  for  all  storms.  To  compare  the  results  of 
Burke’s  quadratic  fit  with  the  J77  model,  Tx  was  calculated  for  each  GRACE  latitude 
bin  in  2004  using  Burke’s  quadratic  fit  given  in  equations  (6)  and  (7)  and  using  the 
interpolation/iteration  technique  with  J77  model  tables  described  in  section  III.  Resulting 
values  of  Too  were  then  orbit  averaged.  While  the  J77  table  interpolation/iteration 
technique  is  much  more  computationally  intensive  than  the  quadratic  fit,  it  produces  a 
more  accurate  representation  of  the  true  J77  model  output  because  it  includes  a  maximum 
error  tolerance  described  in  section  III. 

Figure  14  shows  the  orbit-averaged  exospheric  temperature  values  resulting  from 
Burke’s  quadratic  fit  (red)  and  the  J77  interpolation/iteration  technique  (black)  for  all  of 
2004.  Burke’s  fit  exhibits  significantly  less  variation  than  the  J77  interpolation/iteration 
method,  especially  when  temperatures  drop  below  850K  or  rises  above  1000K. 
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Figure  14:  Comparison  of  the  orbit-averaged  exospheric  temperature  resulting 
from  Burke’s  quadratic  fit  to  the  J77  model  (red)  and  the  J77 
interpolation/iteration  method  developed  for  this  thesis  (black)  as  a  function  of 
Julian  Date,  2004  where  the  date  is  counted  from  1  January,  2004. 


Figure  15  shows  the  quadratic  fit  results  (  7’0o,Burfce)»  plotted  as  a  function  of  the 
interpolation/iteration  results  (foo./nterp)-  If  the  methods  produced  equivalent  results,  all 
data  would  fall  along  the  line  with  a  slope  of  one  and  a  y-intercept  of  zero  shown  as  a 
black  dashed  line.  Burke’s  fit  only  closely  matches  the  J77  tables  in  a  narrow  range 
around  850K.  Above  this  value  Burke’s  fit  produces  temperatures  lower  than  the  J77 
tables  and  below  this  value  Burke’s  fit  produces  temperatures  significantly  higher  than 
J77  tables.  Because  the  quadratic  fit  does  not  closely  match  the  J77  interpolation/iteration 
technique  over  the  whole  range  of  temperatures,  and  departs  significantly  for  the  low 
temperatures  below  850K  that  are  common  near  solar  min  from  2006-2008,  the  J77 
interpolation/iteration  technique  was  used  to  generate  the  observed  GRACE  temperatures 
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used  in  this  thesis.  Too  was  calculated  for  each  GRACE  latitude  bin  and  converted  to  T\a 


using  Equation  (18).  The  temperature  values  were  then  orbit-averaged. 


Figure  15:  Orbit-averaged  exospheric  temperature  from  Burke’s  J77  quadratic  fit 
(Too , Burke)  plotted  as  a  function  of  orbit-averaged  exospheric  temperature  from 
the  J77  interpolation/iteration  method  ( Too,interp )  Is  shown  with  blue  dots.  Data 
shown  is  from  2004.  The  dotted  black  line  has  a  slope  of  1  and  a  y-intercept  of  0. 


Model  Results 

Using  the  observed  orbit-average  T\a  values  calculated  from  GRACE  data  via  the 
J77  interpolation/iteration  technique,  the  model  was  run  using  methods  one,  two  and 
three  described  in  section  III  for  each  storm  period  listed  in  Table  1.  Table  2  shows  the 
model  results.  The  columns  from  left  to  right  list  the  year  of  the  storm,  the  starting  day  of 
the  storm,  the  minimum  value  of  the  Dst  index  during  the  storm  period,  the  value  of  the 
F  io.7a  index  on  the  first  day  of  the  storm,  the  value  of  the  pre-storm  equilibrium 
temperature  (T \a),  and  the  change  in  T i/2uv  over  the  storm  period  (AT i/2uv)  defined  as 
the  difference  between  the  value  of  Ti/2uv  at  the  end  of  the  storm  period  and  Ti/2°.  All  of 
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these  values  are  the  same  for  each  method.  Next  values  for  the  coupling  constant  (a)  with 


units  of 


- ,  relaxation  constant  (x)  with  units  of  hrs,  and  the  relative  T1/2  RMS  error 

hr-mV 


(RMS)  in  percent  resulting  from  each  method  for  each  storm  are  listed. 


Table  2:  Model  Results 


Start 

Min 

rj,  0 

1 1/2 

A 

Method  1 

Method  2 

Method  3 

Year 

Day 

Dst 

Fl0.7a 

T 1/2UY 

a 

T 

RMS 

a 

T 

RMS 

a 

T 

RMS 

2002 

246 

-109 

179.0 

1214.4 

-0.07 

39.47 

6.80 

1.087% 

39.61 

6.78 

1.088% 

39.61 

6.79 

1.087% 

2002 

250 

-181 

179.1 

1273.2 

21.58 

67.28 

3.64 

1.828% 

67.82 

3.55 

1.884% 

68.38 

3.32 

2.332% 

2002 

272 

-176 

175.7 

1134.5 

-5.49 

37.93 

5.95 

1.932% 

37.11 

6.12 

1.921% 

36.40 

6.18 

1.901% 

2002 

276 

-146 

175.0 

1141.2 

31.96 

28.65 

7.64 

0.745% 

27.22 

7.43 

0.750% 

31.33 

4.97 

1.018% 

2002 

296 

-98 

174.4 

1158.4 

18.04 

31.95 

4.23 

1.980% 

38.61 

3.40 

2.066% 

36.54 

3.55 

2.619% 

2002 

324 

-128 

174.0 

1116.8 

-24.16 

25.44 

4.25 

2.631% 

26.79 

4.18 

2.654% 

27.02 

5.30 

2.040% 

2003 

149 

-144 

128.8 

983.3 

-49.18 

41.42 

3.20 

3.715% 

37.96 

3.68 

3.563% 

36.02 

4.39 

2.174% 

2003 

229 

-148 

127.7 

937.2 

-10.75 

41.35 

6.98 

3.599% 

42.59 

6.85 

3.626% 

41.79 

7.43 

3.444% 

2003 

324 

-422 

129.2 

983.4 

17.93 

71.14 

4.41 

3.939% 

69.68 

4.45 

3.960% 

70.56 

4.18 

4.486% 

2004 

22 

-149 

126.0 

1044.2 

-32.74 

26.29 

5.08 

2.122% 

26.73 

5.26 

2.084% 

26.29 

6.46 

1.763% 

2004 

94 

-112 

120.8 

885.1 

-22.38 

37.23 

7.49 

1.452% 

36.87 

7.71 

1.403% 

35.46 

8.23 

1.468% 

2004 

204 

-197 

106.4 

876.0 

-117.63 

31.45 

7.16 

4.456% 

30.72 

7.82 

4.279% 

35.32 

10.03 

3.003% 

2004 

208 

-197 

107.7 

856.2 

-26.60 

33.99 

6.72 

1.603% 

34.47 

6.84 

1.586% 

33.62 

7.47 

1.554% 

2004 

243 

-126 

109.1 

793.8 

7.17 

44.63 

5.58 

2.270% 

43.32 

5.72 

2.331% 

41.13 

5.76 

2.338% 

2004 

312 

-373 

109.9 

970.0 

-10.76 

61.12 

7.05 

1.739% 

61.07 

7.11 

1.746% 

60.05 

7.34 

1.691% 

2004 

314 

-289 

110.3 

1122.4 

-84.39 

63.56 

2.96 

2.731% 

59.20 

3.34 

2.746% 

44.58 

5.22 

1.941% 

2005 

127 

-127 

93.5 

914.9 

28.97 

46.38 

5.88 

2.354% 

48.74 

5.47 

2.476% 

50.14 

5.11 

3.048% 

2005 

135 

-263 

94.0 

904.3 

3.49 

52.88 

7.33 

2.070% 

52.73 

7.28 

2.056% 

52.81 

7.19 

2.077% 

2005 

148 

-138 

94.2 

793.7 

6.43 

26.89 

11.96 

2.242% 

27.39 

11.63 

2.153% 

29.48 

10.77 

2.164% 

2005 

163 

-106 

95.3 

830.2 

-37.58 

30.27 

8.27 

3.027% 

30.13 

9.04 

2.989% 

28.84 

11.08 

2.610% 

2005 

236 

-216 

94.8 

802.7 

-31.30 

46.59 

5.33 

3.401% 

47.46 

5.37 

3.430% 

44.87 

6.72 

1.926% 

2005 

243 

-131 

94.7 

797.1 

-19.41 

47.31 

5.33 

2.669% 

46.68 

5.55 

2.586% 

45.55 

6.08 

2.076% 

2006 

98 

-80 

81.2 

839.3 

-8.75 

22.99 

11.51 

1.548% 

24.20 

11.11 

1.508% 

24.19 

12.09 

1.483% 

2006 

103 

-111 

81.5 

816.6 

-12.25 

31.61 

9.09 

2.587% 

31.62 

9.20 

2.543% 

31.31 

9.78 

2.336% 

2006 

348 

-146 

80.0 

788.2 

-19.73 

35.52 

8.81 

2.104% 

35.35 

9.18 

2.103% 

34.44 

9.84 

1.961% 

2007 

142 

-63 

75.8 

725.7 

-16.93 

20.98 

8.11 

2.312% 

21.88 

7.92 

2.336% 

20.32 

10.41 

1.686% 

2007 

191 

-39 

74.6 

758.6 

0.63 

22.51 

12.28 

0.830% 

22.43 

12.28 

0.825% 

22.43 

12.18 

0.826% 

2007 

195 

-45 

74.5 

767.6 

-7.36 

29.67 

8.74 

1.058% 

30.39 

8.71 

1.041% 

29.99 

9.35 

0.973% 

2007 

218 

-34 

74.0 

709.3 

-3.11 

25.63 

9.16 

0.973% 

25.58 

9.26 

1.005% 

24.87 

9.83 

0.940% 

2007 

298 

-52 

71.5 

715.8 

-1.67 

30.55 

5.58 

1.571% 

30.69 

5.56 

1.564% 

30.90 

5.59 

1.521% 

2007 

323 

-63 

69.9 

697.7 

-1.70 

24.09 

12.96 

1.973% 

23.89 

13.15 

1.967% 

23.79 

13.22 

1.953% 

2007 

351 

-38 

70.5 

821.4 

-30.61 

22.20 

5.25 

2.470% 

22.28 

5.48 

2.467% 

24.32 

7.58 

1.425% 

2008 

31 

-44 

70.3 

708.2 

-5.25 

18.16 

13.73 

1.308% 

18.40 

13.64 

1.253% 

18.64 

13.99 

1.181% 

2008 

68 

-72 

70.4 

711.4 

1.24 

17.99 

11.46 

1.391% 

18.53 

11.05 

1.404% 

18.56 

11.21 

1.423% 

2008 

86 

-43 

70.8 

731.7 

-13.79 

23.28 

14.90 

0.964% 

24.74 

14.19 

0.905% 

25.89 

14.53 

0.791% 

2008 

166 

-40 

70.5 

714.7 

-3.87 

28.66 

8.46 

0.776% 

28.81 

8.49 

0.798% 

27.80 

9.29 

0.693% 

2008 

194 

-40 

69.9 

681.9 

-1.89 

25.57 

7.42 

0.793% 

25.79 

7.46 

0.792% 

25.70 

7.70 

0.783% 

2008 

247 

-51 

69.3 

645.1 

-3.59 

18.47 

12.39 

1.204% 

18.24 

12.73 

1.198% 

17.83 

13.28 

1.161% 
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Table  3  lists  the  mean  and  standard  deviation  of  a,  x,  and  the  relative  Tm  RMS 


error  for  each  method.  The  mean  values  of  a,  t  and  relative  RMS  error  are  very  similar 


between  methods  one  and  two.  The  only  difference  between  these  two  methods  is  the 
factor  in  Equation  (32).  Since  is  calculated  using  the  24-hr  change  in  the 

F10.7  index  via  Equations  (33)  and  (34)  it  is  generally  a  small  term  resulting  in  little 
difference  between  methods  one  and  two.  Method  three  has  mean  alpha  values  lower  than 
methods  one  and  two  and  mean  tau  values  slightly  higher.  The  mean  relative  RMS  error 


for  method  three  is  10%  less  than  that  of  method  one. 


Table  3:  Model  Statistics 


Mean 

a 

Std  Dev 

a 

Mean 

T 

Std  Dev 

T 

Mean 

RMS 

Std  Dev 

RMS 

Lowest  R 

Number 

MS  Storms 

Percentage 

Method  1 

35.03 

13.90 

7.71 

3.07 

2.04% 

0.94% 

7 

18% 

Method  2 

35.15 

13.50 

7.74 

3.01 

2.03% 

0.93% 

4 

11% 

Method  3 

34.65 

12.87 

8.25 

3.05 

1.84% 

0.80% 

27 

71% 

Figure  16  shows  histograms  of  relative  T\n  RMS  error  values  (top),  a  values 
(bottom  left),  and  x  values  (bottom  right)  for  each  of  the  three  methods.  Relative  RMS 
error  values  for  methods  one  and  two  have  a  diffuse  peak  from  1.5  -  2.5%,  and  a  range  of 
0.76%  -  4.46%.  The  standard  deviation  is  very  similar  for  the  two  methods:  0.94%  for 
method  one  and  0.93%  for  method  two.  Method  three  has  a  stronger  relative  RMS  error 
peak  between  1.5%  -  2%  but  a  larger  overall  range  from  0.69%  -  4.49%.  In  general  the 
method  three  errors  are  more  tightly  packed,  with  a  standard  deviation  of  0.80%. 
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Figure  16:  Histograms  showing  relative  RMS  errors  (top),  a  values  (bottom 
left),  and  x  values  (bottom  right)  for  method  one  (blue),  method  2  (green),  and 
method  3  (maroon).  Bars  for  methods  1-3  indicate  values  which  fall  between 
adjacent  labels  on  the  x  axis.  For  example,  the  bottom  left  histogram  shows  that 
methods  1,  2,  and  3  each  had  3  storms  with  a  values  between  10  and  20. 


Coupling  constant  values  range  from  roughly  17-68  for  all  methods  with  values 
most  frequently  falling  between  20  and  30.  Higher  a  values  amplify  the  impact  of  8vs  on 
T 1/2  in  the  model  due  to  the  term  asvs(t)  in  Equations  (31),  (32),  and  (38)  for  methods 
one,  two,  and  three  respectively.  Therefore,  storms  with  higher  temperature  rises  will 
require  higher  alpha  values  in  order  to  model  them  accurately,  a  values  for  method  three 
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are  slightly  more  closely  spaced  than  those  of  methods  one  and  two  as  evidenced  by  the 
slightly  smaller  standard  deviation  of  12.87  for  method  three  vs.  13.90  for  method  one 
and  13.50  for  method  two. 

The  relaxation  constant  controls  how  quickly  T \a  recovers  to  near  pre-storm 
levels  after  sVs  returns  to  normal.  Storms  with  a  faster  recovery  result  in  lower  x  values. 
Relaxation  constant  values  ranged  from  roughly  3  -  1 5  for  all  methods  with  values  falling 
most  frequently  between  5  -  7.5.  The  range  in  x  values  was  significantly  smaller  than  the 
range  in  a  values  indicating  that  the  recovery  period  of  geomagnetic  storms  is  less 
variable  than  the  main  phase.  The  spread  in  x  values  was  similar  for  all  methods,  with 
standard  deviations  just  over  3  hours. 

It  should  be  noted  that  the  seemingly  small  difference  in  mean  T  \a  errors  between 
method  three  (1.84%)  and  method  one  (2.04%)  is  significant  due  to  its  impact  on  density 
errors.  Small  changes  (or  errors)  in  thermospheric  temperatures  result  in  large  changes 
(or  errors)  in  thermospheric  densities.  A  brief  example  from  the  J77  model  illustrates  this 
point.  If  the  observed  exospheric  temperature  is  taken  to  be  700K,  a  1.84%  error  in  Too 
(matching  the  mean  method  three  error)  would  generate  a  density  error  of  13.58%  at 
an  altitude  of  500  km.  If  instead  the  2.04%  Ti/2  error  from  method  one  were  applied,  the 
resulting  density  error  would  increase  to  15.14%  at  500  km.  In  this  case  the  0.2% 
increase  in  temperature  error  produces  a  1.56%  increase  in  density  error  illustrating  that 
the  slight  increase  in  temperature  accuracy  produced  by  method  three  is  operationally 
relevant.  Density  results  will  be  discussed  in  further  detail  later  in  the  document. 
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Individual  Storms 


Figure  17  shows  the  T\a  results  of  the  model  for  the  CME  storm  from  Julian  Date 
(JD)  204-210,  2004.  This  storm  period  actually  includes  three  distinct  CMEs  hitting  the 
earth  in  rapid  succession  as  evidenced  by  the  magnetospheric  electric  field  data  shown  in 
the  bottom  plot.  The  start  of  the  storm  period  is  defined  as  the  time  the  first  CME  hits  on 
JD  204,  indicated  by  the  vertical  red  line.  The  second  and  third  CMEs  can  be  seen  in  the 
abrupt  rises  in  electric  field  magnitude  on  JD  206  and  just  prior  to  JD  209.  For  this  storm, 
method  three  was  significantly  better  than  methods  one  or  two,  producing  a  relative 
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Figure  17:  Model  results  for  the  CME  storm  from  Julian  Date  204-210,  2004.  The 
top  plot  shows  observed  GRACE  T1/2  (red  dots),  along  with  model  T1/2  results  for 
methods  one  (black),  two  (pink),  and  three  (green).  The  dotted  red  line  shows  the 
pre-storm  equilibrium  temperature  for  methods  one  and  two  and  the  black  dotted 
line  shows  the  results  of  the  approximation  Ty2  ~  T \/2uv  used  for  method  3.  The 
bottom  plot  shows  the  electric  field  values  calculated  from  ACE  data  as  a  function  of 
time.  The  red  vertical  line  indicates  the  storm  starting  time. 
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T 1/2  RMS  error  of  3.00%  compared  with  4.46%  for  method  one  and  4.28%  for  method 
two.  Method  three  clearly  outperformed  the  other  two  methods  especially  in  fitting  the 
T 1/2  peak  from  the  first  CME  right  around  JD  205  and  during  the  T\a  minimum  just  prior 
to  the  third  CME  at  the  end  of  JD  208. 

The  large  differences  between  method  three  and  methods  one  and  two  for  this 
storm  are  due  to  the  relatively  large  change  in  T1j2uv  during  the  storm  period  from  a  pre¬ 
storm  equilibrium  value  of  876.00  K  down  to  758.37  K  by  the  end  of  the  storm  as  shown 
by  the  dotted  black  line  in  Figure  17.  Method  one  ignores  this  change  entirely  by  treating 
T1/2uv  as  a  constant  value  through  the  storm.  Method  two  takes  the  change  into  account 

partially  through  the  term  in  equation  (32)  but  does  not  allow  Tj°/2  to  vary  with 

T1/2uv  during  the  storm  period.  For  this  storm,  while  the  overall  AT1/2uv  is  -1 17.63K 
over  the  storm  period,  the  rate  of  change  remains  small,  never  dropping  below 

-1.79  K/hr.  By  not  allowing  Ty2  to  vary  with  T1/2uv  during  the  storm,  methods  one  and 
two  result  in  an  awkward  situation  on  JD  208,  when  observed  GRACE  Tl/2  values  drop 
below  the  supposed  UV  contribution  to  7l/2-  This  means  that  if  methods  one  and  two 
were  to  be  accurate  during  this  period,  they  would  have  to  produce  a  negative  value  for 
T1/2sw>  the  amount  of  temperature  rise  due  to  the  solar  wind,  which  is  an  unphysical 
result. 

Method  three  avoids  this  situation  by  taking  into  account  the  change  in  T1/2lJV 
during  the  storm  period  by  approximating  Ty2  ~  7lj2UV  which  results  in  Equation  (38). 
Allowing  T1/2UV  to  decrease  through  the  storm  period  by  definition  (Equation  (26)) 
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increases  T^sw-  This  results  in  a  higher  a  value  for  method  three  for  this  storm,  35.32, 
than  methods  one  and  two,  31.45  and  30.72,  respectively.  Similarly,  the  decrease  in 
T1/2uv  during  the  storm  period  results  in  an  increased  x  value  for  method  three,  10.03, 
compared  with  methods  one  and  two,  7.16  and  7.82,  respectively.  The  decreasing  Tlj2\jv 
plays  a  role  similar  to  the  relaxation  constant  and  helps  the  modeled  T i/2  recover  after  svs 
decreases.  Since  the  decreasing  T1/2UV  performs  a  similar  role  to  the  relaxation  constant, 
method  three  results  in  a  higher  x  value. 

Figure  18  shows  the  results  of  methods  one,  two  and  three  for  the  CIR  storm  from 
JD  351-356,  2007.  The  overall  Ty2  increase  for  this  storm  over  the  pre-storm  equilibrium 
value  of  82 IK  is  about  60K.  Similarly  to  the  CME  storm  in  Figure  17,  this  CIR  has  a 
decreasing  T i/2Uv  throughout  the  storm  period.  Methods  one  and  two  produce  virtually 

identical  results,  due  to  the  fact  that  the  for  this  storm  is  very  small,  never 

dropping  below  -0.32  K/hr.  Method  three  accounts  for  the  drop  in  Ti/2uv  of  about  30K 
resulting  in  a  relative  RMS  error  of  only  1.43%  compared  with  2.47%  for  methods  one 
and  two.  The  drop  in  Ty2Uv  also  results  in  higher  a  and  x  values  compared  with  methods 
one  and  two  for  the  same  reasons  as  the  CME  storm. 

Method  three  produced  larger  errors  than  method  one  or  method  two  for  1 1 
storms  in  the  sample.  Figure  19  shows  the  results  for  one  of  these  storms,  the  CME-storm 
from  JD  250-252,  2002.  For  this  storm  method  1  produced  the  lowest  RMS  error.  Unlike 
the  two  previous  storms,  in  this  CME  T i/2uv  increases  throughout  the  storm  period.  This 
causes  x  for  method  three  to  be  lower  than  methods  one  or  two,  resulting  in  the 
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Figure  18:  Model  results  for  the  CIR  storm  from  JD  351-356,  2007.  The  top  plot 
shows  observed  GRACE  Tm  (red  dots),  along  with  model  Ti/2  results  for  methods  one 
(black),  two  (pink),  and  three  (green).  The  dotted  red  line  shows  the  pre-storm 
equilibrium  Ti/2  for  methods  one  and  two  and  the  black  dotted  line  shows  the  results 
of  the  approximation  Ty2  «  T1j2uv  used  for  method  3.  The  bottom  plot  shows 
magnetospheric  electric  field  values.  The  red  vertical  line  indicates  the  storm  start 
time. 


decreased  accuracy  near  the  peak  of  the  storm.  In  addition,  method  three  models  the 
recovery  phase  of  the  storm  worse  than  methods  one  or  two  because  the  increasing  Ti/2uv 

(XT 

increases  the  method  three  model  1/2  (t)  in  Equation  (37)  during  a  time  when  observed 
T 1/2  is  decreasing  . 
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Figure  19:  Model  results  for  the  CME  from  JD  250-252,  2002.  The  top  plot  shows 
observed  GRACE  T1/2  (red  dots),  along  with  model  T2/2  results  for  methods  one 
(black),  two  (pink),  and  three  (green).  The  dotted  red  line  shows  the  pre-storm 
equilibrium  T1/2  for  methods  one  and  two  and  the  black  dotted  line  shows  the  results 
of  the  approximation  T\j2  ~  T \/2uv  used  for  method  three.  The  bottom  plot  shows 
magnetospheric  electric  field  values.  The  red  vertical  line  shows  the  storm  start  time. 


Method  Comparison 

Overall,  method  three  produced  the  lowest  errors  most  frequently  among  the  38 
storms  tested.  Table  3  shows  that  method  three  produced  the  lowest  relative  RMS  error 
for  27  of  the  38  storms  studied  (71%),  while  methods  one  and  two  only  had  the  lowest 
error  for  7  storms  (18%)  and  4  storms  (11%),  respectively.  The  method  that  produced  the 
lowest  error  for  a  given  storm  was  strongly  dependent  on  the  nature  of  the  change  (either 
increasing  or  decreasing)  in  Ti/2uv  over  the  storm  period. 

Figure  20  shows  the  model  method  that  produced  the  lowest  relative  Ti/2  RMS 
error  as  a  function  of  the  change  in  T i/2uv  (AT  i/2uv)  over  the  storm  period.  AT  i/2uv  was 
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Figure  20:  Model  method  producing  the  lowest  relative  Ti/2  RMS  error  as  a 
function  of  the  change  in  T1/2uv  ( AT|  2r\  )  over  the  storm  period. 


defined  as  the  difference  between  the  value  of  Ti/2uv  at  the  end  of  the  storm  period  and 
the  pre-storm  equilibrium  temperature.  The  majority  of  the  storms  (28  of  38)  had  a 
decrease  in  T i/2uv  over  the  storm  period.  Method  three  produced  the  lowest  error  for  27 
of  those  28  storms.  A  negative  ATi/2Uv  allows  method  three  to  accurately  characterize  the 
recovery  period  of  the  storm  with  a  larger  x  value  than  methods  one  and  two  because 
some  of  the  decrease  in  Ty2  is  accounted  for  by  the  decreasing  Ti/2Uv  at  the  end  of  the 
storm  period.  This  larger  x  value  in  turn  means  that  the  peak  of  the  storm  is  more 

accurately  modeled  because  a  larger  x  tends  to  increase  ~~  (t).  Equation  (37),  during 
the  growth  phase  of  the  storm  when  T i/2uv  is  still  near  the  pre-storm  equilibrium  level. 
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The  single  storm  with  a  negative  AT  i/2uv  that  was  not  modeled  best  by  method  three  was 
the  CIR  storm  from  JD  94-98,  2004.  For  this  storm  Ti/2uv  rose  above  the  pre-storm 
equilibrium  temperature  on  JD  95-96  before  dropping  below  on  JD  97.  Method  two  was 
the  best  method  for  this  storm. 

Method  three  did  not  produce  the  lowest  error  for  any  of  the  10  storms  with  a 
positive  ATi/2uv.  Storms  with  higher  values  of  ATi/2uv  tended  to  be  modeled  by  method 
one  best  while  all  three  of  the  storms  with  positive  ATi/2uv  values  for  which  method  two 
produced  the  best  results  had  ATi/2uv<  6.5K.  Positive  ATi/2uv  forced  the  x  value  for 
method  three  to  decrease,  resulting  in  an  underestimate  of  the  peak  T1/2  values  of  the 
storm.  In  addition,  the  increase  in  Ti/2uv  caused  method  three  to  model  the  recovery 
phase  of  the  storm  worse  than  method  one.  Both  of  these  problems  are  clearly  illustrated 
in  the  JD  250,  2002  CME  shown  in  Figure  19. 

Comparison  with  Burke,  2011 

Burke,  2011  determined  a  values  for  37  of  the  38  storms  listed  in  Table  1.  To 
compare  Burke’s  2011  results  with  current  results,  Burke’s  a  values  need  to  be  divided 
by  a  storm-average  value  of  the  conversion  factor  D,  from  Equation  (19),  to  account  for 
the  fact  that  he  modeled  Too  instead  of  T1/2.  The  storm  average  value  of  D  was  generally 

,  yy\ 

near  0.95.  After  conversion,  Burke’s  (2011)  results  have  a  mean  a  value  of  36.07 - 

v  '  hr-mV 

and  a  standard  deviation  of  17  Km  .  Burke  (201 1)  treated  Toouv  as  a  constant  in  his 

hr-mV 

model  similarly  to  method  one  here.  Burke’s  mean  a  and  standard  deviation  of  a  are 
higher  than  those  resulting  from  method  one,  shown  in  Table  3. 
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There  are  several  reasons  for  the  difference.  First,  Burke  (201 1)  assumed  x  was 
constant  over  all  storms.  The  value  of  6.84  hrs  (after  conversion)  was  determined  from 
the  single  storm  of  JD  204-210,  2004  by  Burke  et  al.,  2009  using  a  different  calibration  of 
GRACE  data.  Allowing  x  to  change  between  storms  impacts  the  value  of  a.  In  addition, 
Burke,  2011  used  one -hour  averaged  ACE  data  to  calculate  svs  values  and  a  time-step  of 
1  hr  when  applying  Equation  (31)  rather  than  the  1  minute  time  step  used  here.  Further, 
Burke  2011  used  slightly  different  start  and  end  times  for  the  storm  periods  than  used 
here,  used  the  quadratic  fit  to  the  J77  model  in  Equations  (6)  and  (7)  to  calculate 
observed  To0;  and  used  a  different  method  to  calculate  the  pre-storm  equilibrium 
temperature;  namely  using  the  value  of  Too  at  the  start  time  of  the  storm  rather  than 
averaging  over  the  12  hours  prior.  Finally,  Burke,  201 1  used  trial  and  error  to  determine 
the  best  value  for  a  in  an  effort  to  align  modeled  Too  with  the  observed  peak.  In  this 
thesis,  the  Nelder-Mead  simplex  direct  search  method  was  applied  in  order  to  rigorously 
determine  the  best  a  and  x  values.  Out  of  all  these  procedural  differences,  the  method  of 
determining  a  and  x  values  for  each  storm  has  the  biggest  impact  on  results. 

To  quantify  the  impact  of  using  a  rigorous  method  to  determine  optimum  a  and  x 
values  the  model  was  re-run  for  the  storm  of  JD  204-210,  2004  using  Equation  (31)  from 
method  one  and  matching  Burke’s  procedures  as  closely  as  possible.  One -hour  ACE  data 
was  used  to  calculate  svs  and  a  one-hour  time  step  was  used  in  Equation  (31).  In  addition 
Too  was  modeled  instead  of  T 1/2,  Burke’s  pre-storm  equilibrium  temperature  was  used 
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(Too°=  810K),  and  Burke’s  quadratic  fit  and  orbit-averaging  technique  was  used  to 
determine  observed  GRACE  Too. 

Figure  21  shows  a  contour  plot  of  relative  T»  RMS  errors  as  a  function  of  a  and  x 
values  for  the  JD  204-210,  2004  storm.  Contour  plots  were  generated  by  running  the 
model  (method  one  here)  for  a  100x100  grid  of  a  and  x  values  and  computing  the  relative 
RMS  error  resulting  from  each  a  and  x  pair.  The  a  and  x  values  from  Burke,  2011  (point 
A)  result  in  a  relative  Tx  RMS  error  of  3.82%  compared  with  a  relative  Too  RMS  error  of 
2.58%  resulting  from  optimal  a  and  x  values  (point  C)  determined  by  using  the 


a  [K*m*hr'1*mV1] 

Figure  21:  Contour  plot  of  relative  RMS  errors  (%)  in  resulting  from  different  a  and  x 
values  using  the  procedures  from  Burke,  2011  for  the  CME  storm  from  JD  204-210,  2004. 
Point  A  shows  the  a  and  x  values  reported  by  Burke,  2011.  Point  B  shows  the  a  and  x 
values  that  result  from  using  the  Nelder-Mead  simplex  direct  search  method  to  minimize 
the  error  between  the  peak  observed  GRACE  T*,  and  the  modeled  value.  Point  C  shows 
the  a  and  x  values  resulting  from  using  the  Nelder-Mead  simplex  direct  search  method  to 
minimize  the  relative  RMS  error  in  over  the  entire  storm  period. 
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Nelder-Mead  simplex  direct  search  method  to  minimize  the  relative  RMS  error.  Point  B 
shows  the  a  and  x  values  resulting  from  the  Nelder-Mead  method  applied  to  minimize 
the  error  between  the  peak  observed  GRACE  Too  and  the  modeled  value.  This  is  the  error 
Burke  (2011)  was  trying  to  minimize  via  trial  and  error. 

The  contour  plot  shows  that  the  relative  RMS  error  is  a  relatively  shallow 

function  of  x  within  the  range  of  5  -  8  hrs  and  a  within  the  range  of  30  to  45 

"TTZ 

Because  of  this  the  difference  between  Burke’s  a  and  x  values  (44.00 - ,6.50  hrs)  and 

v  hr-mV’  J 

"T7X 

the  optimal  values  (31.83  ,7.25  hrs)  only  decreases  relative  RMS  error  from  3.82% 

to  2.58%.  This  shows  that  while  Burke’s  trial  and  error  method  came  close  to  the  optimal 
values,  the  more  rigorous  approach  produces  superior  results  and  helps  explain  the 
difference  between  the  method  one  a  values  and  Burke’s  results.  The  fact  that  relative 
RMS  error  is  a  relatively  shallow  function  of  a  and  x  near  the  minimum  also  suggests 
that  it  is  possible  to  produce  acceptable  results  with  a  and  x  values  different  than  the 
optimal  values. 

Solar  Cycle  Dependence  of  a  and  x 

Burke,  2011  suggested  that  the  values  for  a  might  vary  throughout  the  solar  cycle 
as  a  function  of  Fio.7a.  Since  the  relative  RMS  error  is  a  shallow  function  of  a  and  x  near 
the  minimum,  it  is  reasonable  to  expect  that  the  model  will  produce  low  errors  with 
values  of  a  and  x  that  are  different  than  the  optimal  values  for  each  storm.  In  order  to  test 
this,  best- fit  curves  were  constructed  to  produce  a  and  x  as  functions  of  the  Fio.7a  value, 
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the  162-day  running  average  of  the  F10.7  index  used  in  the  J77  model,  at  the  start  of  each 
storm  period  shown  in  Table  2.  Since  FI 0.7a  is  a  162-day  average,  it  changes  very  little 
over  the  1-6  day  storm  periods  used  in  this  Thesis.  Least-Squares  fits  were  constructed 
using  data  from  method  three  because  it  proved  to  produce  the  lowest  errors  of  the  three 
methods  for  most  of  the  storms  in  the  sample. 

Figure  22  shows  method  three  a  values  as  a  function  of  Fioja  for  all  storms.  A 
linear  fit,  shown  in  black,  produces  a  tenuous  correlation  of  R  =  0.21.  Robinson  and 
Yondrak,  1984  showed  that  both  the  ion-electron  production  rate,  which  impacts  particle 
precipitation,  and  ionospheric  conductance,  which  impacts  joule  heating,  depend  on 
yjF-i o.7  (Burke,  2011).  Since  a  accounts  for  the  energy  transfer  from  the  magnetosphere 
to  the  thermosphere  via  joule  heating  and  particle  precipitation  in  the  driven-dissipative 
model,  it  is  reasonable  to  fit  a  as  a  quadratic  function  of  ^  F107a  (Burke,  2011).  The 
quadratic  fit  of  a  to  A/F10  7a  ,  shown  in  red,  produces  a  correlation  of  R  =  0.40  which  is 
much  improved  over  the  linear  fit.  Figure  23  shows  x  as  a  function  of  Fi0.7a  for  all  storms 
using  method  three.  Again,  a  quadratic  least-squares  fit  was  constructed  of  x  as  a  function 


of 


.  The  correlation  of  R 


0.53  is  higher  than  the  correlation  for  a. 
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Figure  22:  Coupling  constant  (a)  as  a  function  of  Fi0.7a  shown  with  blue  dots.  The  black 
line  and  text  show  the  best  linear  fit  to  the  data.  The  red  line  and  text  show  the  best  fit 
of  a  as  a  quadratic  function  of  J F10  7a . 


Figure  23:  Relaxation  constant  (x)  as  a  function  of  Fi0.7a  shown  with  blue  dots.  The 
black  line  shows  the  best  fit  of  x  as  a  quadratic  function  of  Jf10  7cl. 
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Outliers 


Figure  22  highlights  the  fact  that  three  storms  have  significantly  higher  values  of 
a  than  the  others.  The  JD  250,  2002  CME  has  an  a  value  of  68.38,  the  JD  324,  2003 
CME  has  an  a  value  of  70.56  and  the  JD  3 12,  2004  CME  has  an  a  value  of  60.05.  All 
other  storms  have  a  <  53.  The  high  a  values  for  the  2003  and  2004  storms  are  due  to  the 
fact  that  they  were  by  far  the  strongest  storm  in  the  sample.  For  the  2003  storm,  the 
minimum  Dst  value  was  an  extreme  -422  nT  and  it  had  Ty2  rise  of  536  K  over  the  pre¬ 
storm  equilibrium  value.  The  2004  storm  had  a  minimum  Dst  value  of  -389  nT  and  a  Ti/2 
rise  of  538K.  These  Ti/2  rises  are  270K  larger  than  the  next  highest  rise  in  the  sample. 
The  high  a  values  for  the  JD  324,  2003  CME  and  the  JD  3 14,  2004  CME  are  a  result  of 
the  large  storm-time  rise  in  Ti/2. 

In  contrast,  the  JD  250,  2002  storm  has  a  minimum  Dst  of  -181  nT  and  a  storm¬ 
time  rise  in  T i/2  of  21  OK.  While  this  is  a  large  rise  in  T i/2,  there  were  1 1  storms  in  the 
sample  which  had  a  larger  T i/2  rise  yet  a  smaller  a.  The  JD  250,  2002  CME  had  such  a 
high  a  value  because  its  storm  time  T  y2  rise  resulted  from  a  relatively  weak  svs  signature 
with  a  maximum  of  1.37  mV/m  shown  in  Figure  19.  For  comparison,  the  JD  204,  2004 
storm  shown  in  Figure  17  had  a  T  i/2  rise  18%  higher  than  JD  250,  2002  yet  the  maximum 
s vs  value  was  58%  higher.  The  JD  250,  2002  storm  had  a  temperature  rise  that  was 
disproportionately  larger  than  the  solar  wind  energy  contribution,  modeled  with  the  asvs 
term,  would  indicate.  For  this  storm,  the  extra  energy  came  from  a  spike  in  solar  EUV 
energy  shown  in  Figure  24. 
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Figure  24  plots  five-minute  average  EUV  (26-34  nm)  flux  (black  line),  as 
measured  by  the  Solar  and  Heliospheric  Observatory  (SOHO)  satellite,  as  a  function  of 
JD,  2002  for  the  time  period  of  the  JD  250,  2002  CME.  Daily  F10.7  (solid  red  line  and  x’s) 
and  F  10.7a  (dotted  red  line  with  x’s)  values  are  also  shown.  The  blue  vertical  line  indicates 
the  storm  start  time.  A  large  spike  in  EUV  flux,  due  to  a  solar  flare,  is  clearly  seen  just 
after  the  storm  start  time.  Since  the  model  accounts  for  EUV  energy  with  the  daily  F10.7 
index,  it  does  not  capture  variations  on  such  short  time  scales  as  the  flare  seen  during  this 
storm.  As  a  result,  the  model  has  to  account  for  this  EUV  flare  energy  by  attributing  it  to 
the  solar  wind  contribution,  asvs-  Since  £vs  is  small  for  this  storm,  the  only  way  to 
increase  the  solar  wind  contribution  is  by  increasing  the  a  value. 


Julian  Date,  2002 

Figure  24:  Five-minute  average  Solar  EUV  flux  (26-34  nm)  measured  by  the  SOHO 
satellite  (black)  shown  as  a  function  of  JD,  2002.  Daily  F10.7  (solid  red  line  and  x’s)  and 
F10.7a  values  (dotted  red  line  and  x’s)  are  also  shown.  The  blue  line  indicates  the  start 
time  of  the  JD  250,  2002  CME  storm. 
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Impact  of  Storm  Type 

Since  CME-driven  storms  have  features  distinctly  different  than  CIR-driven 
storms,  their  response  to  solar  cycle  changes  may  not  be  equivalent.  To  test  this,  least- 
squares  fits  of  a  and  x  to  ^/F107a  were  constructed  for  each  storm  type  separately.  Of  the 
38  storms  in  the  sample,  25  were  CIR  storms  and  13  were  CME  storms.  Table  1  lists  the 
storm  type  of  each  storm.  Figure  25  shows  a  as  a  quadratic  function  of  F107a  for  CME 
storms  (blue)  and  CIR  storms  (red).  The  value  of  a  for  CME  storms  exhibits  a  very  low 
correlation,  R  =  0.12,  and  is  nearly  a  straight  line.  CME  storms  generally  had  higher  a 
values  than  CIR  storms.  In  fact  all  storms  with  a  >  50  are  CMEs  while  all  of  the  storms 
with  a  <  25  are  CIR  storms.  CIR  storms  are  fit  much  better  as  a  function  of  ■yjF1  a7a  with 
a  correlation  of  R  =  0.60,  much  improved  from  the  all-storms  fit.  The  CIR  storms  are  fit 
better  because  both  CIR  occurrence  and  F10.7  index  both  exhibit  a  27-day  period  linked  to 
solar  rotations.  Since  CME  occurrence  is  irregular  with  no  characteristic  spacing,  CMEs 
are  not  as  well  correlated  with  F10.7  measurements  (Borovsky  and  Denton,  2006). 

Figure  26  shows  x  as  a  quadratic  function  of  A/F10  7a  for  CME  storms  (blue)  and 
CIR  storms  (red).  Both  best  fit  curves  are  similar,  with  correlations  of  R=0.47  and 
R=0.5 1  for  CIRs  and  CMEs,  respectively.  The  correlations  are  worse  for  both  storm 
types  than  the  correlation  of  the  all  storms  fit  indicating  that  x  is  not  strongly  dependent 
on  storm  type.  In  general,  x  values  are  higher  for  CIR  storms.  All  storms  with  x  >  1 1  are 
CIR-driven. 
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Figure  25:  Coupling  constant  (a)  as  a  function  of  Fi0.7a.  Best  fits  of  a  as  a  quadratic 
function  of  J F10  7a  are  shown.  CME  storms  are  shown  in  blue,  CIR  storms  in  red. 


Figure  26:  Relaxation  constant  (x)  as  a  function  of  F10.7a.  Best  fits  of  x  as  a  quadratic 
function  of  ^jF10  7a  are  shown.  CME  storms  are  shown  in  blue,  CIR  storms  in  red. 
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Best  Fit  a  and  x  Results 


All  38  storms  were  run  with  all  storm  fit  a  and  x  values  determined  with  the 
quadratic  equations  shown  in  Figure  22  and  Figure  23  using  method  three  procedures.  In 
addition,  CME  storms  were  run  using  the  CME  fit  a  and  x  values  and  CIR  storms  were 
run  using  CIR  fit  a  and  x  values  determined  using  functions  shown  in  Figure  25  and 
Figure  26.  Table  4  shows  the  mean  Ty2  relative  RMS  errors,  calculated  using  Equation 
(39),  that  result  from  method  three  using  optimal  a  and  x  values  for  each  storm  listed  in 
Table  2,  along  with  relative  RMS  errors  that  result  from  the  all  storms  fit,  CME  fit,  and 
CIR  fit  a  and  x  values.  Column  2,  labeled  All  Storms  Mean,  shows  that  the  mean  Ti/2 
RMS  error  increased  from  1.84%  to  3.15%  for  all  38  storms  when  using  the  all  storms  fit 
a  and  x  values.  For  CME  and  CIR  storms,  errors  increased  from  the  optimal  values  when 
both  the  all  storms  fit  and  the  storm  specific  fit  was  applied.  For  both  storm  types,  the 
storm  type  specific  fit  values  of  a  and  x  produced  a  lower  average  RMS  error  than  the  all 
storms  fit.  Applying  best-fit  curves  to  determine  a  and  x  for  each  storm  also  created  more 
spread  in  the  relative  Ti/2  RMS  error  values. 

Table  4:  Relative  Ti/2  RMS  Error  Results  using  Best  Fit  a  and  x  values  with 


Method  Three 


All  Storms 

CME  Storms 

CIR  Storms 

Mean 

Std  Dev 

Mean 

Std  Dev 

Mean 

Std  Dev 

Optimal  a  and  x 

1.84% 

2.17% 

0.83% 

1.67% 

All  Storms  Fit  a  and  x 

3.15% 

1.68% 

4.37% 

2.00% 

2.52% 

CME-Fit  a  and  x 

4.01% 

1.77% 

CIR-Fit  a  and  x 

2.24% 

■KIIM 
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Table  4  shows  standard  deviations  of  the  relative  RMS  error  for  each  of  the 


different  fits.  In  all  cases,  the  standard  deviation  increased  over  the  value  for  the  optimal 
a  and  x  case.  Figure  27  shows  histograms  of  the  relative  T\a  RMS  error  for  all  storms 
(top),  CME  storms  (bottom  left)  and  CIR  storms  (bottom  right).  In  the  all  storm 
histogram  we  see  the  all  storm  fit  error  values  are  spread  over  a  much  larger  range  than 


1  5  2  2.5  3  3.5  4  4.5 


Relative  RMS  Error  (%) 


Figure  27:  Histograms  showing  relative  Ti/2  RMS  errors  from  method  three  using 
best-fit  a  and  x  values  for  all  storms  (top),  CME  storms  (bottom  left),  and  CIR 
storms  (bottom  right).  Errors  for  optimal  a  and  x  values  are  shown  in  green,  errors 
for  all  storm  fit  a  and  x  values  are  shown  in  gray,  errors  for  CME-fit  a  and  x 
values  are  shown  in  blue,  and  errors  for  CIR-fit  a  and  x  values  are  shown  in  red. 
Bars  show  the  number  of  storms  which  resulted  in  a  relative  Ti/2  RMS  error 
between  the  adjacent  values  of  the  x-axis. 
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the  optimal  error  values.  The  maximum  error  using  optimal  a  and  x  values  lies  between 
4  -  4.5%  while  the  maximum  error  using  all  storm  fit  a  and  x  values  lies  between 
8  -  8.5%.  Similar  patterns  are  seen  for  CME  storms  and  CIR  storms.  The  spread  is  most 
significant  for  CME  storms,  where  the  standard  deviation  more  than  doubles  from  0.83% 
for  the  optimal  case  to  2.00%  for  the  all  storm  fit  case. 

Individual  Storms 

Table  5  shows  the  results  of  the  all  storms  fit  and  the  CME  fit  a  and  x  values 
applied  to  the  CME  storm  of  JD  204-210,  2004.  The  optimal  a  and  x  for  this  storm  are 
included  for  comparison.  In  this  case,  the  CME  fit  produced  a  higher  error  than  the  all¬ 
storms  fit.  This  is  not  surprising  as  the  correlation  for  the  CME-fit  function  for  a  was 
very  low  (R  =  0.12).  a  and  x  values  resulting  from  the  all  storms  fit  and  the  CME  fit  are 
very  similar  for  this  storm,  resulting  in  the  similar  model  T 1/2  curves  for  these  two  cases 
seen  in  Figure  28  as  the  solid  pink  (all  storms  fit)  and  black  (CME  fit)  lines.  Both  the 
best  fit  a  and  x  values  produce  model  T1/2  curves  that  are  below  the  optimal  one  (shown 
in  green)  resulting  in  decreased  accuracy  especially  during  the  second  and  third  T1/2 
peaks  on  JD  207  and  JD  209.  The  best  fit  a  values  are  higher  than  the  optimal  value, 
suggesting  that  model  T 1/2  increase  more  rapidly  when  sys  increases,  however  the  best  fit 
x  values  are  lower  than  the  optimal  case  which  indicates  a  quicker  recovery  time  and 
decreases  the  modeled  increase  in  T1/2.  For  this  storm,  the  decrease  in  x  wins  out  and 
causes  the  model  T 1/2  for  the  all  storm  and  CME  fit  cases  to  lag  below  the  optimal  case 
during  the  peaks  on  JD  205  and  207. 
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Table  5: 

Results  1 

for  JD  204-210,  2004  CME 

a 

T 

Relative  T|/2  RMS  Error 

Optimal 

35.32 

10.03 

3.00% 

All  Storms  Fit 

41.93 

7.74 

3.30% 

CME  Fit 

42.08 

7.17 

3.76% 

203  204  205  206  207  208  209  210 

Julian  Date,  2004 


Figure  28:  Model  results  for  the  CME  storm  from  Julian  Date  204-210,  2004.  The  top 
plot  shows  observed  GRACE  T1/2  (red  dots),  along  with  method  three  model  T1/2  using 
optimal  a  and  x  values  (green),  all  storms  fit  a  and  x  values  (pink),  and  CME  fit  a  and 
x  values  (black).  The  the  black  dotted  line  shows  Ti/2UV.  The  bottom  plot  shows  the 
electric  field  values  calculated  from  ACE  data  as  a  function  of  time.  The  red  vertical 
line  indicates  the  storm  starting  time. 


A  contour  plot  of  relative  T\n  RMS  errors  (%)  as  a  function  of  a  and  x  is  shown 
in  Figure  29  for  the  JD  204-210,  2004  CME  storm.  As  expected,  relative  RMS  error  is  a 
relatively  shallow  function  of  a  and  x  around  the  minimum,  shown  as  point  A 
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corresponding  to  the  optimal  a  and  x  values  for  this  storm.  Points  B  and  C  correspond  to 
a  and  x  values  calculated  using  the  all  storms  fit  and  the  CME  fit,  respectively.  Even 
though  these  values  are  not  very  close  to  the  minimum,  the  relative  RMS  error  for  both  is 
still  less  than  four  percent. 


Figure  29:  Contour  plot  of  relative  T1/2  RMS  errors  (%)  as  a  function  of  a  and  t  for  the  JD 
204-210,  2004  CME  storm  using  method  three.  Point  A  shows  the  optimal  a  and  t  values, 
point  B  shows  the  all  storms  fit  a  and  t  values,  and  point  C  shows  the  CME  fit  a  and  t 
values. 


Table  6  shows  the  results  of  the  all  storms  fit  and  the  CIR  fit  a  and  x  values 
applied  to  the  CIR  storm  of  JD  351-356,  2007.  The  optimal  values  of  a  and  x  are 
included  for  comparison.  For  this  storm  the  CIR- fit  produced  a  lower  relative  RMS  error 


than  the  all-storms  fit.  Best  fit  a  values  for  this  storm  are  very  close  to  the  optimal  value 
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while  best  fit  x  values  are  3  -  3.5  hrs  higher  than  the  optimal  value.  The  increased  x  value 


indicates  that  the  best  fit  models  should  result  in  higher  Ti /2  values,  especially  after  sys 
decreases,  because  higher  x  values  equate  to  a  longer  e-fold  recovery  time. 


Table  6:  Results  for  JD  351-356, 2007  CIR 


a 

X 

Relative  Ti/2  RMS  Error 

Optimal 

24.32 

7.58 

1.43% 

All  Storms  Fit 

23.84 

11.26 

2.76% 

CIR  Fit 

23.32 

10.53 

2.24% 

Figure  30  shows  the  model  results  for  the  JD  351-356,  2007  CIR  storm  and  as 
expected,  both  the  all  storm  fit  and  the  CIR  fit  a  and  x  values  have  T  m  curves  that  are 
higher  than  the  optimal  case  resulting  in  the  increased  errors  shown  in  Table  6.  As  with 
the  JD  204,  2004  CME,  here  the  all  storm  fit  and  the  storm  specific  fit  a  and  x  values  are 
similar  leading  to  the  small  difference  between  the  all  storm  and  CIR  fit  T\a  curves. 

A  contour  plot  of  relative  T\a  RMS  error  (%)  as  a  function  of  a  and  x  is  shown 
for  the  JD  351-356,  2007  CIR  in  Figure  31.  Point  A  shows  the  location  of  the  minimum 
error  resulting  from  optimum  a  and  x  values  while  points  B  and  C  show  the  locations  of 
the  all  storm  fit  and  CIR  fit  a  and  x  values,  respectively.  Again,  the  error  is  a  relatively 
shallow  function  of  a  and  x  around  the  minimum.  Any  a,  x  pair  within  the  ranges  of  20  < 
a  <  28  and  5  <  x  <  10  results  in  a  relative  RMS  error  of  less  than  3%  for  this  storm. 
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Figure  30:  Model  results  for  the  CIR  storm  from  Julian  Date  351-356,  2007.  The  top 
plot  shows  observed  GRACE  T2/2  (red  dots),  along  with  method  three  model  Ti/2  using 
optimal  a  and  x  values  (green),  all  storms  fit  a  and  x  values  (pink),  and  CIR  fit  a  and 
x  values  (black).  The  black  dotted  line  shows  Ti/2UV.  The  bottom  plot  shows  the 
electric  field  values  calculated  from  ACE  data  as  a  function  of  time.  The  red  vertical 
line  indicates  the  storm  starting  time. 


General  Applicability 

The  results  of  the  model  using  best  fit  a  and  x  values  indicate  that  relatively  low 
errors  can  be  obtained  using  model  parameters  determined  without  prior  knowledge  of 
storm-time  T1/2  values.  To  test  this  conjecture  the  model  with  best-fit  a  and  x  values  was 
applied  to  two  storms,  one  CME  and  one  CIR,  outside  the  original  sample  of  38  storms 
that  were  used  to  determine  the  best  fit  a  and  x  curves.  Due  to  the  constraints  of  GRACE 
data  availability,  both  test  storms  were  within  the  same  time  frame  of  the  original  storms 
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Figure  31:  Contour  plot  of  relative  Tj/2  RMS  errors  (%)  as  a  function  of  a  and  x  for  the  JD 
351-356,  2007  CIR  storm  using  method  three.  Point  A  shows  the  optimal  a  and  x  values,  point 
B  shows  the  all  storms  fit  a  and  x  values,  and  point  C  shows  the  CME  fit  a  and  x  values. 


(2002-2008).  Test  storms  were  selected  based  on  the  availability  of  ACE  data  needed  to 
calculate  storm  time  svs  values.  Table  7  shows  from  left  to  right  the  year,  start  time,  end 
time,  minimum  Dst  index  during  the  storm  period,  Fio.7a  value  on  day  one  of  the  storm, 
the  pre-storm  equilibrium  temperature  (T  m),  AT  i/2uv  for  each  storm  period,  and  the 


storm  type  for  each  of  the  test  storms. 


Table  7:  Test  Storm  Data 


Year 

Storm  Start 

Day  Hour  Min  Sec 

Storm  End 

Day  Hour  Min  Sec 

Min 

Dst 

Fl0.7a 

Tl/2° 

A  T1/2UV 

Storm 

Type 

2003 

308  3  36  0 

309  6  0  0 

-69 

132.06 

1049.4 

-24.39 

CME 

2004 

42  4  48  0 

45  0  0  0 

-109 

123.46 

928.5 

-20.16 

CIR 

89 


Table  8  shows  the  results  for  the  JD  308-309,  2003  CME.  Although  the  relative 


T 1/2  RMS  error  for  both  the  all  storms  fit  and  CME-fit  a  and  x  values  are  more  than 
double  the  optimal  error,  they  are  also  less  than  the  average  error  for  CME  storms  in  the 
original  sample.  This  indicates  that  the  best  fit  a  and  x  values  are  reasonable  even  outside 
the  original  sample. 


Table  8: 

Results  1 

for  JD  30! 

8-309,  2003  CME 

a 

X 

Relative  T|/2  RMS  Error 

Optimal 

39.01 

4.46 

0.98% 

All  Storms  Fit 

45.54 

6.47 

2.89% 

CME  Fit 

46.42 

5.50 

2.17% 

Figure  32  shows  the  T1/2  curves  resulting  from  the  optimal,  all  storms  fit,  and 
CME  fit  a  and  x  values.  Both  fits  cause  the  model  to  overestimate  the  peak  T1/2  value  and 
do  not  recover  as  fast  as  the  observed  T1/2.  The  CME  fit  has  slightly  lower  errors  than  the 
all  storms  fit  due  to  its  lower  x  value,  which  causes  T1/2  to  drop  faster  during  the  recovery 
period  and  close  the  gap  between  the  modeled  and  observed  T1/2. 

Figure  33  shows  a  contour  plot  of  the  relative  Ti/2  RMS  error  as  a  function 
of  a  and  x  for  the  JD  308,  2003  storm.  The  errors  for  this  storm  are  more  sensitive  to 
changes  in  x  compared  to  the  storms  shown  previously.  However,  within  the  range  of  3  < 
x  <  6,  a  can  take  any  value  between  35  and  48  and  still  produce  an  error  of  less  than  3%. 
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Figure  32:  Model  results  for  the  CME  storm  from  Julian  Date  308-309,  2003.  The  top 
plot  shows  observed  GRACE  T1/2  (red  dots),  along  with  method  three  model  Ti/2  using 
optimal  a  and  x  values  (green),  all  storms  fit  a  and  x  values  (pink),  and  CME  fit  a  and 
x  values  (black).  The  black  dotted  line  shows  Ti/2uv*  The  bottom  plot  shows  the 
electric  field  values  as  a  function  of  time.  The  red  vertical  line  shows  storm  start  time. 
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Figure  33:  Contour  plot  of  relative  Ti/2  RMS  errors  (%)  as  a  function  of  a  and  x  for  the 
JD  308-309,  2003  CME  storm  using  method  three.  Point  A  shows  the  optimal  a  and  x 
values,  point  B  shows  the  all  storms  fit  a  and  x  values,  and  point  C  shows  the  CME  fit  a 
and  x  values. 
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Results  from  the  CIR  test  storm  are  shown  in  Table  9.  Again,  the  relative  Ty2 


RMS  error  for  both  the  all  storms  fit  and  the  CIR  fit  a  and  x  values  are  less  than  or  equal 


to  the  average  errors  for  CIR  storms  in  the  original  sample.  The  CIR  storm  fit  produces 
lower  errors  than  the  all  storms  fit  for  this  storm.  The  T y2  curves  produced  by  each  model 
run  are  shown  in  Figure  34.  The  CIR  fit  model  actually  matches  the  peak  Ti/2  value  just 


Table  9:  Results  for  JD  42-45,  2004  CIR 


a 

T 

Relative  T1/2  RMS  Error 

Optimal 

28.74 

8.56 

1.05% 

All  Storms  Fit 

44.99 

6.81 

2.52% 

CIR  Fit 

42.46 

5.98 

1.43% 

41.5  42  42.5  43  43.5  44  44.5  45 

Julian  Date,  2004 


Figure  34:  Model  results  for  the  CIR  from  JD  42-45,  2004.  The  top  plot  shows 
observed  T2/2  (red  dots),  along  with  method  three  model  T1/2  using  optimal  a  and  x 
values  (green),  all  storms  fit  a  and  x  values  (pink),  and  CIR  fit  a  and  x  values  (black). 
The  black  dotted  line  shows  T1/2UV.  The  bottom  plot  shows  the  electric  field  values 
calculated  from  ACE  data.  The  red  vertical  line  indicates  the  storm  start  time. 
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after  JD  42.5  more  accurately  than  the  optimal  case.  The  increase  in  error  is  due  to  the 
modeled  Tin  over-reacting  to  the  second  svs  peak  on  JD  43  resulting  in  over-estimates  of 
T 1/2.  The  higher  a  and  x  values  of  the  all  storms  fit  cause  the  all-storm  fit  model  to 
produce  higher  T  m  at  all  times  compared  with  the  CIR  fit.  A  notable  feature  of  this 
storm  is  the  fact  that  the  best  fit  a  and  x  values  are  significantly  different  than  the  optimal 
values  yet  the  relative  RMS  errors  do  not  increase  drastically.  Figure  35  shows  the 
relative  T1/2  RMS  error  as  a  function  of  a  and  x.  There  is  a  very  broad  range  of  a  and  x 
values  which  result  in  errors  of  less  than  3%  for  this  storm  and  both  best  fit  models  fall 
within  the  range. 


a  value 

Figure  35:  Contour  plot  of  relative  Ti/2  RMS  errors  (%)  as  a  function  of  a  and  x  for  the  JD  42- 
45,  2004  CIR  storm  using  method  three.  Point  A  shows  the  optimal  a  and  x  values,  point  B 
shows  the  all  storms  fit  a  and  x  values,  and  point  C  shows  the  CIR  fit  a  and  x  values. 
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Density  Errors 

In  order  to  compare  the  results  of  the  driven-dissipative  model  to  published 
HASDM  density  errors,  model  T i/2  values  are  used  to  calculate  model  densities.  Using 
the  methods  described  in  section  III,  model  orbit-averaged  densities  were  computed  for 
each  storm  in  the  38-storm  sample  and  relative  density  RMS  errors  for  each  storm  were 
calculated.  In  general,  higher  temperature  errors  should  result  in  higher  density  errors. 
Figure  36  shows  the  relative  RMS  error  in  model  orbit-average  density  plotted  as  a 
function  of  relative  RMS  error  in  model  orbit-average  Ti/2.  In  general,  the  errors  follow 
the  expected  trend  with  high  Ti/2  errors  resulting  in  high  density  errors.  However,  there 
are  three  outlier  storms  with  density  errors  greater  than  13%  resulting  from  Ti/2  errors  of 
less  than  1.5%. 


Figure  36:  Model  relative  orbit-average  density  RMS  error  plotted  as  a  function  of 
relative  orbit-average  Ti/2  RMS  errors.  Blue  dots  show  storms  from  2002  -  JD  290, 
2007.  Red  x’s  show  storms  from  JD  290,  2007  -  2008. 
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The  outlier  storms  are  caused  in  part  by  differences  in  the  orbit-average 
techniques  used  when  computing  observed  GRACE  orbit-averaged  T\a  and  model  orbit- 
averaged  density.  Figure  37  illustrates  the  two  different  techniques  as  applied  to 
exospheric  temperature.  Observed  GRACE  orbit-averaged  T\a  was  calculated  via  the 
bin-averaging  technique  where  the  J77  model  is  applied  to  calculate  a  temperature  in 
each  latitude  bin  prior  to  orbit-averaging.  Producing  observed  orbit-averaged  temperature 
values  via  the  bin-  averaging  technique  is  mathematically  preferable  to  the  whole-orbit 
technique  because  the  latter  technique  is  akin  to  calculating  an  average  of  averages. 

These  techniques  are  not  equivalent  and  produce  slightly  different  results.  The  average 
relative  RMS  difference  between  the  orbit-averaged  T,_  produced  by  the  two  techniques 
was  small  for  most  storms,  ranging  between  0.20%  and  3.07%  with  an  average  of 
1.33%.  Figure  38  shows  the  results  of  the  two  techniques  for  the  storm  with  the  largest 
difference  between  the  two.  Temperature  values  calculated  via  the  bin-averaging 
technique  (blue)  are  lower  than  those  generated  from  whole-orbit  technique  (green)  for 
this  storm  and  all  storms  in  the  sample. 


Bin-AveragingTechnique 


Whole-OrbitTechnique 


Figure  37:  Diagram  outlining  two  different  orbit-averaging  techniques,  p  is  the  mass  density,  H  is  the 
height  above  sea  level,  and  T,  is  the  exospheric  temperature.  Orbit-average  values  are  indicated  by 
p,  H,  and  Tm. 
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Figure  38:  Orbit-Average  exospheric  temperature  as  a  function  of  Julian  Date,  2008. 

Too  values  calculated  via  bin-averaging  are  shown  in  blue.  T„o  values  calculated  via  the 
whole-orbit  technique  are  shown  in  green. 

Ideally,  to  ensure  consistency,  the  bin-averaging  technique  would  be  applied  to 
calculated  model  orbit-averaged  densities  from  modeled  Tm  values.  This  would  require 
calculation  of  a  model  density  in  each  latitude  bin  and  then  the  results  would  be  averaged 
over  entire  orbits.  As  discussed  in  section  III,  the  current  model  is  not  formulated  to 
accurately  produce  exospheric  temperatures  in  specific  latitude  bins.  Therefore,  it  is  not 
meaningful  to  apply  the  bin-averaging  technique  to  produce  orbit-averaged  model 
densities.  Instead,  the  whole-orbit  technique  is  used  to  calculate  model  orbit-averaged 
densities  from  model  orbit-averaged  Tx  and  observed  orbit-averaged  height  values. 

The  result  of  this  mixing  of  techniques  is  a  built-in  model  density  error  caused  by 
the  difference  between  the  two  techniques.  If  the  model  performed  perfectly  and 
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produced  model  Tm  values  equivalent  to  the  observed  T\a  values  calculated  using  the 
bin-averaging  technique,  the  orbit-average  model  densities  resulting  from  the  whole-orbit 
technique  would  be  lower  than  the  observed  orbit-averaged  densities.  Since  the  difference 
between  the  Too  produced  by  the  two  techniques  is  small  for  most  storms,  the  resulting 
density  errors  are  not  contaminated  significantly.  However,  for  two  of  the  outlier  storms 
(JD  86,  2008  and  JD  247,  2008)  the  relative  Too  RMS  difference  between  the  two 
techniques  was  greater  than  2.4%  resulting  in  model  density  values  that  were 
significantly  lower  than  the  observed  values  despite  the  fact  that  the  T 1/2  errors  for  these 
storms  were  quite  small. 

The  third  outlier  (JD  351,  2007)  had  a  relatively  small  difference  between 
exospheric  temperatures  calculated  with  the  two  different  orbit  averaging  techniques 
which  suggests  that  there  is  another  factor  influencing  the  density  errors  for  outlier 
storms.  All  of  the  outlier  storms  fell  at  the  end  of  the  sample  period,  at  the  end  of  2007 
and  into  2008.  This  matches  the  time  period  of  the  last  solar  minimum  (solar  cycle 
23/24),  which  was  centered  in  November,  2008  (Emmert,  et  al.,  2010).  EUV  flux  and 
thermospheric  density  during  the  last  solar  min  were  markedly  lower  than  all  five  other 
solar  mins  observed  since  the  start  of  the  space  age  (Solomon  et  al.,  2010).  This  is  likely 
to  impact  the  driven-dissipative  model  because  it  is  based  on  the  J77  model,  which  was 
built  using  fits  to  observed  data  from  satellite  drag  measurements.  In  addition,  the  driven- 
dissipative  model  accounts  for  EUV  flux  by  using  the  F10.7  index.  During  the  last  solar 
min,  observed  EUV  flux  decreased  by  15%  compared  with  averages  from  the  previous 
five  solar  mins  while  F10.7  values  were  only  down  by  5%  (Chen  et  al.,  201 1).  This 


97 


indicates  a  change  in  the  relationship  between  F10.7  and  EUV  flux  which  could  impact  the 
model  through  Equation  (33). 

Emmert,  et  al.,  2010,  studied  whether  the  changes  in  thermospheric  densities 
during  the  last  solar  min  could  be  modeled  by  solely  changing  exospheric  temperatures. 

A  model  density  profile  was  constructed  using  a  Bates-Walker  diffusive  equilibrium 
profile  (Walker,  1965).  This  type  of  profile  is  similar  to  the  one  used  in  the  J77  model. 
Perturbing  exospheric  temperature  alone  did  not  result  in  a  model  density  profile  which 
matched  the  average  density  profile  observed  in  2008  -  2009.  Instead,  changes  to 
exospheric  temperature  and  thermospheric  composition  were  both  needed  to  produce  a 
model  density  profile  which  matched  observations  (Emmert,  et  al.,  2010).  This  suggests 
that  the  driven-dissipative  model,  which  only  perturbs  temperatures,  will  not  perform 
well  during  the  last  solar  min. 

Emmert  et  al.,  2010  found  that  the  difference  between  observed  thermospheric 
density  departed  from  1986  -  2007  climatology  by  more  than  10%  beginning  in 
November,  2007.  To  eliminate  the  impact  of  the  last  solar  min  on  model  density  results, 
all  storms  from  November,  2007  through  2008  (shown  as  red  x’s  in  Figure  36)  were 
discarded.  This  removed  all  outlier  storms  and  resulted  in  the  expected  trend  of 
increasing  T 1/2  errors  resulting  in  increasing  density  errors  as  shown  by  the  blue  dots  in 
Figure  36. 

After  discarding  storms  during  the  last  solar  min  density  errors  resulting  from  the 
driven-dissipative  model  compare  favorably  with  HASDM.  This  is  significant  because 
the  model  results  include  the  built-in  error  introduced  by  the  difference  between  orbit- 
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averaging  techniques.  Table  10  shows  the  mean  and  standard  deviation  of  the  relative 
density  RMS  errors  for  the  29  storms  from  the  original  sample,  defined  in  Table  1, 
occurring  from  2002  -  October,  2007.  Results  are  shown  for  all  storms  and  separated  by 
storm  type  for  the  optimal  a  and  x  values  as  well  as  the  all  storms  fit  a  and  x  values  and 
the  storm  type  (CME  or  CIR)  fit  a  and  x  values.  The  optimal  and  all  storms  fit  mean 
density  errors  of  11.18%  and  18%,  respectively,  compare  favorably  with  the  mean 
HASDM  storm-time  error  of  17%.  Optimal  errors  for  the  CME  and  CIR  storms  are  also 
well  below  HASDM  errors.  Density  errors  resulting  from  best  fit  a  and  x  values  for  CIR 
storms  are  less  than  HASDM  mean  errors  while  for  CMEs,  best  fit  a  and  x  values  result 
in  mean  density  errors  slightly  higher  than  mean  HASDM  storm-time  errors. 


Table  10:  Relative  Density  RMS  Error  resulting  from  Best  Fit  a  and  x  Values  for 

29  storms  from  2002  -  October  2007 


CIR  Storms 

Mean 

Std  Dev 

Optimal  a  and  x 

11.18% 

3.52% 

11.69% 

2.72% 

10.77% 

4.10% 

All  Storms  Fit  a  and  x 

18.00% 

9.78% 

22.28% 

11.96% 

14.53% 

5.96% 

CME-Fit  a  and  x 

20.23% 

7.51% 

CIR-Fit  a  and  x 

12.47% 

5.23% 

Figure  36,  along  with  comparisons  between  mean  T\a  errors  shown  in  Table  4 
and  mean  density  errors  shown  in  Table  10,  illustrates  that  for  any  given  storm  density 
errors  are  much  higher  than  temperature  errors.  The  increase  in  error  is  due  to  the  fact 
that  small  changes  in  temperature  result  in  large  changes  in  density  within  the  J77  model. 
Figure  39  shows  the  observed  orbit-averaged  density  for  the  JD  204,  2004  CME  (red 
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dots)  along  with  model  densities  resulting  from  optimal  a  and  x  values  (green  line),  all 
storm  fit  a  and  x  values  (blue  line),  and  CME  fit  a  and  x  values  (black  line).  Similar  to 
the  T1/2  results  shown  for  this  storm  in  Figure  28,  the  all  storm  fit  and  CME  fit  a  and  x 
values  result  in  densities  that  are  lower  than  both  observed  and  optimal  model  values  for 
the  second  and  third  peaks  on  JD  207  and  209.  For  this  storm  the  observed  orbit-average 
densities  range  from  2.17  x  10~16  to  10.72  x  10~16  g/cm3,  an  increase  of  394%.  In  contrast, 
observed  orbit-average  T i/2  values  for  this  storm,  shown  in  Figure  28,  range  from  841.18 
to  1 124.4K,  an  increase  of  only  34%.  In  other  words  relatively  small  changes  (or  errors) 
in  temperature  values  result  in  large  changes  (or  errors)  in  density  values.  This  sensitivity 
explains  the  difference  between  model  temperature  and  model  density  errors. 


Julian  Date,  2004 


Figure  39:  Model  density  results  for  the  CME  storm  from  Julian  Date  204-210,  2004. 
Observed  GRACE  orbit-average  density  is  shown  by  red  dots,  along  with  model 
density  values  resulting  from  optimal  a  and  x  values  (green  line),  all  storm  fit  a  and  x 
values  (blue  line),  and  CME  fit  a  and  x  values  (black  line). 
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Model  density  errors  were  also  computed  for  the  two  test  storms  outside  of  the 
original  sample,  defined  in  Table  7.  Table  1 1  shows  the  relative  density  RMS  error 
resulting  from  optimal,  all  storm  fit,  and  storm  type  (CME  or  CIR)  fit  a  and  x  values  for 
these  two  storms.  As  expected,  density  errors  are  higher  than  Tm  errors  but  still  low 
compared  to  the  mean  HASDM  storm-time  density  error  of  17%. 


Table  11:  Relative  Density  RMS  Error  resulting  from  Best  Fit  a  and  x  Values  for 

Two  Test  Storms 


JD  308-309,  2003 

JD  42-45,  2004 

CME 

CIR 

Optimal 

4.51% 

6.20% 

All  Storms  Fit 

11.86% 

15.38% 

Storm  Type  Fit 

8.19% 

8.97% 
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V.  Conclusions  and  Recommendations 


Conclusions  of  Research 

This  project  has  produced  conclusions  in  three  main  areas.  First,  the  impact  of  the 
EUV  term  in  the  driven-dissipative  model  was  explored.  It  was  found  that  over  the  entire 
38  storm  sample,  method  three  procedures  (allowing  both  the  time  rate  of  change  of 

Ti/2uv,  d7’1^w,  and  approximating  Ti/2uv(t)  =  T 1/2°)  were  on  average  slightly  more 

accurate  than  methods  one  and  two.  Method  three  had  a  mean  relative  Ti/2  RMS  error 
1.84%  compared  to  2.04%  and  2.03%  for  methods  one  and  two,  respectively.  The  impact 
of  the  treatment  of  the  UV  contribution  was  strongly  dependent  on  the  character  of  the 
change  in  Ti/2uv  over  the  storm  period.  Method  three  produced  the  smallest  relative  T1/2 
RMS  error  for  27  of  the  28  storms  in  the  sample  that  had  decreasing  Ti/2uv  profiles.  In 
contrast  method  one,  which  treated  Ty2uv  as  a  constant,  produced  the  smallest  relative 
T1/2  RMS  error  for  seven  of  the  10  storms  with  increasing  T1/2UV  profiles.  Method  two, 

which  allowed  to  vary  but  treated  T1/20  as  a  constant,  produced  results  very 

similar  to  method  one  for  all  storms.  In  general,  for  the  declining  phase  of  the  solar  cycle, 
the  full  treatment  of  the  UV  contribution  used  in  method  three  is  the  most  accurate 
variation  of  the  driven-dissipative  model. 

The  second  conclusion  of  this  thesis  is  that  the  two  empirical  parameters  in  the 
driven-dissipative  model  exhibit  solar  cycle  dependence  and  can  be  determined  as 
quadratic  functions  of  yjF10Ja,  where  the  Fio.7a  value  used  is  the  value  at  the  start  of  the 
storm  period.  This  is  important  because  it  provides  a  way  to  determine  model  parameters 
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without  any  prior  knowledge  of  the  storm  to  be  modeled.  Using  best-fit  model  parameters 
degraded  model  accuracy  slightly,  but  the  mean  relative  T\n  RMS  error  for  all  storms 
remained  small  at  less  than  3.2%.  Around  the  minimum,  relative  T\a  RMS  error  is  a 
shallow  function  of  the  model  parameters  allowing  departures  from  the  optimal  values 
without  significantly  increasing  errors.  Model  accuracy  was  improved  slightly  by 
separating  storms  by  type  (CME  or  CIR)  and  determining  model  parameters  separately  as 
functions  of  y/F1 0.7a  for  each  storm  type.  The  general  applicability  of  the  model  and  the 
model  parameter  fits  was  established  by  applying  them  to  two  test  storms  outside  the 
original  sample  of  38  which  resulted  in  errors  similar  to  those  within  the  original  sample. 

The  final  conclusion  of  this  research  is  that  the  driven-dissipative  model,  as 
formulated  in  method  three,  can  be  used  in  conjunction  with  the  J77  model  to  produce 
model  density  values  with  accuracies  similar  to  those  currently  produced  by  HASDM. 
Mean  relative  density  RMS  errors  for  the  model  averaged  11.18%  when  optimal  model 
parameters  were  used  and  18%  when  the  model  parameters  determined  by  the  all-storms 
fit  functions  of  were  used.  These  values  compares  favorably  to  HASDM’ s  mean 

error  of  17%  during  geomagnetic  storming  conditions.  Of  course  the  current  formulation 
of  the  driven-dissipative  model  uses  observed  solar  wind  data  as  the  driver,  which  helps 
produce  accurate  results.  Still,  the  comparison  with  HASDM  is  significant  because  while 
the  driven-dissipative  model  does  not  use  observed  density  data  to  correct  the  model  in 
near  real-time  as  HASDM  does,  it  is  still  able  to  produce  comparable  density  errors.  This 
suggests  that  future  research  could  further  improve  the  driven-dissipative  model. 
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Recommendations  for  Future  Research 


In  addition  to  producing  several  satisfying  conclusions,  this  project  suggests 
several  promising  avenues  of  future  research.  First,  the  set  of  38  storms  used  in  this  thesis 
all  occurred  during  the  declining  phase  of  the  solar  cycle  from  2002-2008.  Since  the 
results  indicate  that  the  most  accurate  method  of  treating  Ti/2uv  depends  strongly  on  the 
character  of  the  Ti/2uv  change  over  the  storm  period  it  would  be  useful  to  expand  the 
storm  sample  to  cover  an  entire  solar  cycle.  This  would  likely  result  in  a  storm  sample 
that  is  more  evenly  split  between  storms  with  increasing  and  decreasing  Ti/2Uv  changes 
and  provide  a  more  rigorous  test  of  the  three  methods  of  treating  T  i/2uv  used  in  this  thesis. 
Currently,  the  ability  to  test  storms  over  an  entire  solar  cycle  is  limited  by  the  availability 
of  the  GRACE  data  used  as  ground  truth  in  the  model. 

A  second  avenue  of  research  is  related  to  the  current  formulation’s  use  of  the  J77 
model  as  a  basis.  The  J77  model  was  used  in  the  current  formulation  to  be  consistent  with 
Burke’s  earlier  work  (2009,  2011).  Current  cutting  edge  models  such  as  FLASDM  and 
JB08  are  based  on  Jacchia’s  1970  model  instead  of  J77.  This  suggests  that  it  may  be 
useful  to  reformulate  the  driven-dissipative  model  to  use  J70  as  a  basis  instead  of  J77. 
Doing  this  should  provide  two  main  advantages.  First,  it  would  allow  experimentation 
with  different  formulas  to  account  for  the  UY  contribution  to  thermospheric  temperature. 
Using  the  J70  model  as  a  basis  would  allow  the  driven-dissipative  model  to  easily  use  the 
J70,  HASDM,  or  JB08  formulations  for  Touv  which  could  result  in  improved  accuracy. 
Second,  using  the  J70  model  as  a  basis  could  allow  the  method  or  the  output  from  the 
driven-dissipative  model  to  be  integrated  with  HASDM  and/or  JB08  in  order  move  from 
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a  model  temperature  to  a  model  density  value.  Doing  so  could  eliminate  the  built-in  error 
that  results  from  the  current  model’s  method  of  converting  from  model  T\n  to  model 
density  and  also  take  advantage  of  the  use  of  observed  data  in  both  HASDM  and  JB08  to 
improve  accuracy. 

A  final  recommendation  for  future  research  would  be  to  move  toward  replacing 
the  observed  solar  wind  data  used  as  a  driver  in  the  current  formulation  with  input  from  a 
current  solar  wind  model.  While  this  would  almost  certainly  degrade  accuracy,  it  is 
necessary  to  make  the  driven-dissipative  model  useful  in  an  operational  sense. 
Determining  how  much  the  use  of  a  model  input  changes  the  results  would  be  important 
in  assessing  the  potential  of  the  driven-dissipative  model  for  use  in  real-world  forecasting 
applications. 
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Appendix  A  -  Solar  Declination  Angle  Calculation 

This  appendix  outlines  the  procedure  from  Meeus,  1991  used  to  calculate  the 
solar  declination  angle.  All  angle  formulas  presented  here  are  in  decimal  degrees.  Let  Y 
denote  the  year  of  the  data  point,  m  the  month  number  of  the  data  point,  and  D  the  day  of 
the  month  including  decimals  (for  example  the  5  th  day  of  the  month  at  12Z  would  mean 
D  =  5.5).  If  in  <2,  replace  Y  with  Y-l  and  m  with  m-2.  Adopting  the  notation  where 
INT( )  denotes  the  integer  part  of  the  argument  within  the  parenthesis,  the  Julian  date 
(JD)  is  calculated  via  the  formula 


JD  =  INT  (365.25(7  +  4716))  +  INT  (30.600  l(m  +  l))  +  D  +  i?-l  524.5  (42) 


Where  B  is  given  by: 


B  =  2- A  +  INT 


(A\ 

UJ 


with 


A  =  INT 


m 

liooj 


Next  a  time  T  is  calculated. 


T  =  JP-2451545  (43) 

36525 

Using  T  the  values  for  the  geometric  mean  longitude  of  the  sun  (L0),  the  mean  anomaly 
of  the  sun  (M),  the  longitude  of  the  ascending  mode  of  the  Moon’s  mean  orbit  on  the 
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ecliptic  (Q),  the  mean  obliquity  of  the  elliptic  (s0),  the  mean  longitude  of  the  sun  (L)  and 


the  mean  longitude  of  the  Moon  (L’)  are  calculated  in  decimal  degrees. 


Lo  =  280.46645  +  36000.7698371  +  0.0003032E2 


(44) 


M  =  357.5291  +  35999.0503J  -  0.0001559J2  -  0.00000048J3 


(45) 


T 3 

Q  =  125. 04452-1934. 13626  ir  +  0.0020708r2  + - 

450000 


(46) 


s 0  =  23 .43929 1 1 1 1  -  0.0 1 3004 1 667F  - 1 .639  x  1 0~7  T2  +  5 .036  x  1 0'7  73  (47) 


L  =  280.4665  +  36000.769877 


(48) 


L'  =  218.3165 +  48 1267. 881377 


(49) 


Using  these  values  a  true  obliquity  (s)  is  calculated. 

s  =  so+As  (50) 


Where  As,  the  nutation  of  the  obliquity,  is  given  by  Equation  (51). 
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(51) 


As  =  2.5556  x  1(T3  cos (Q)  + 1 .5833  x  10~4  cos (2L)  + ... 
2.7778xl0”5  cos(2L')-2.5xlO”5  cos(2fi) 


Next  the  sun’s  equation  of  center  (C)  is  calculated  based  on  the  time  T  and  the  mean 
anomaly  of  the  sun  M 

C  =  (1 .9 1 4600  -  0.0048 17 T-  0.0000 1 4 T2 )  sin(M)  + .. .  (52) 

(0.01 9993  -  0.000 1 0  IT)  sin(2  M)  +  0.000290  sin(3M) 

which  leads  to  the  sun’s  true  longitude  (0). 

®  =  La  +  C  (53) 

Finally,  the  apparent  longitude  of  the  sun,  X,  is  calculated 

/l  =  0-0.00569- 0.00478  sin  (Q)  (54) 

leading  to  the  solar  declination  angle  (8). 

5  =  sin"1  (sin  (e)  sin(/l))  (55) 


108 


Appendix  B  -  The  Nelder-Mead  Simplex  Direct  Search  Method 

MATLAB’s  fminsearch  function  uses  the  Nelder-Mead  simplex  direct  search 
method  (Lagarias,  et  al.,  1998)  to  minimize  a  given  function  by  adjusting  the  specified 
variables.  For  this  research  the  function  to  be  minimized  is  the  RMS  error  function  given 
in  Equation  (39),  denoted  here  as  f(x)  and  the  variables  to  be  adjusted  are  the  coupling 
constant,  a,  and  the  relaxation  constant,  x,  denoted  here  as  the  two  element  vector  x.  The 
number  of  elements  in  the  vector  x  is  denoted  by  n,  here  n=2.  The  algorithm  is  started  by 
providing  initial  values  of  the  vector  x,  denoted  as  x0.  Here  the  initial  values  were  set  at 

"  TTZ 

a°  =  44  —  —  and  r0  =  5.4  hrs.  To  start  the  algorithm  an  initial  simplex  is  created 

around  x0  by  adding  5%  to  each  value  of  x0  one  at  a  time,  resulting  in  n+1  vectors 
defining  the  vertices  of  the  initial  simplex. 

Once  a  simplex  is  defined,  the  vertices  xt  are  ordered  based  on  their  function 
value  such  that  f(xt)  <  /(x2)  <  •••  <  f(xn+1).  During  each  step  in  the  search  iteration 
the  worst  point,  xn+1,  is  discarded  and  replaced  by  a  new  point  via  one  of  the  methods 
outlined  below.  The  iteration  continues  until  the  values  of  the  cost  function  converge  to  a 
user-defined  tolerance.  For  this  research  the  tolerance  was  defined  as  10'6. 

Reflection 

A  reflected  point,  xr,  is  generated  using  the  formula 

xr=2x-xn+l  (56) 

where 
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(57) 


Zn 

X, 

 i- 1  1 


X  = 


/I 


After  the  point  xr  is  generated  the  function  value,  /(xr),  is  calculated.  If 

fOc-i)  <  f(xr)  <  f(xn)  the  point  xr  is  accepted  to  replace  xn+1,  creating  a  new  simplex, 

and  the  iteration  starts  over.  Figure  40  shows  the  simplices  after  a  reflection  step. 


x3 


xr 


Figure  40:  Nelder-Mead  simplices  after  a  reflection.  The  original  simplex  is 
shown  with  a  dashed  line  (Lagarias,  et  al.,  1998). 

Expansion 

If  f(xr )  <  /(xx)  an  expansion  point,  xe>  is  calculated  using  the  formula 

xe=x  +  2{x-xn+l)  (58) 

and  the  resulting  function  value  f(xe)  is  calculated.  If  (xe)  <  f(xr)  ,  the  point  xe  is 
accepted  to  replace  xn+1  and  the  iteration  starts  over.  If  (xe)  >  f(xr)  ,  the  point  xr  is 
accepted  to  replace  xn+1  and  the  iteration  starts  over.  Figure  41  shows  the  simplices  after 
an  expansion  step. 
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Figure  41:  Nelder-Mead  simplices  after  an  expansion.  The  original  simplex  is 
shown  with  a  dashed  line  (Lagarias,  et  al.,  1998). 

Contraction 

If  f(xr)  >  /(xn),  a  contraction  is  performed  between  the  point  x  and  the  point 
xn+1  or  xr  that  produces  the  lowest  function  value.  If  /(xr)  <  /(xn+1)  an  outside 
contraction  is  performed  by  calculating  the  point  xc 


(59) 


2 


and  the  resulting  function  value,  /(xc).  If  /(xc)  <  /(xr),  the  point  xc  is  accepted  to 
replace  xn+1  and  the  iteration  starts  over.  If  /(xc)  >  /(xr),  a  shrink  is  performed  using 
procedures  in  the  next  section.  Figure  42  shows  the  simplices  after  an  outside  contraction. 


Ill 


xr 


xc 


Figure  42:  Nelder-Mead  simplices  after  an  outside  contraction.  The  original 
simplex  is  shown  with  a  dashed  line  (Lagarias,  et  al.,  1998). 


If /(xr)  >  /(xn+1)  an  inside  contraction  is  performed  by  calculating  the  point 

%cc 


Xcc=X  + 


(Xn+l~X) 


(60) 


and  the  resulting  function  value,  /(xcc).  If /(xcc)  <  /(xn+1),  the  point  xcc  is  accepted  to 
replace  xn+1  and  the  iteration  starts  over.  If  /(xcc)  >  /(xn+1),  a  shrink  is  performed 
using  procedures  in  the  next  section.  Figure  43  shows  the  simplices  after  an  inside 
contraction. 


Figure  43:  Nelder-Mead  simplices  after  an  inside  contraction.  The  original 
simplex  is  shown  with  a  dashed  line  (Lagarias,  et  al.,  1998). 
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Shrink 


If  none  of  the  previous  methods  used  to  identify  a  new  point  to  replace  xn+1  in  a 
the  new  simplex  were  successful  a  new  simplex  is  calculated  by  performing  a  shrink 
using  the  formula 


x 


v,-=*i+- 


(61) 


where  i=2...n+l.  The  points  are  ordered  by  increasing  values  of  f{v{)  and  the  new 
simplex  is  defined  by  the  best  point  in  the  old  simplex,  xlf  along  with  the  new  values 
where  again  i  =  2. .  .n+1.  Figure  44  shows  the  simplices  after  a  shrink. 


I 


Xi 


Figure  44:  Nelder-Mead  simplices  after  a  shrink.  The  original  simplex  is  shown 
with  a  dashed  line  (Lagarias,  et  al.,  1998). 


Schematic 

Figure  45  shows  a  schematic  depicting  the  Nelder-Mead  simplex  direct  search 
method  described  above.  The  bold  text  in  each  box  depicts  the  condition  that  must  be 
satisfied  in  order  to  perform  the  action  listed  in  the  box. 
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Figure  45:  Schematic  of  the  Nelder-Mead  simplex  direct  search  method  used  by  MATLAB’s 
fminsearch  function,  described  by  (Lagarias,  et  al.,  1998).  Schematic  is  shown  for  a  function 
of  2  variables,  as  used  in  this  thesis.  The  bold  text  in  each  box  depicts  the  condition  that  must 
be  satisfied  in  order  to  perform  the  action  listed  in  the  box. 
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